基于离心机和数值模拟的软硬互层反倾层状岩质边坡变形特征分析

李彦奇，黄 达，孟秋杰

（河北工业大学土木与交通学院，天津 300400）

在我国西部山区，天然或人工开挖形成了大量的反倾岩质边坡[1−2]。早期工程地质工作者认为反倾层状边坡较为稳定，不易形成贯通面引发滑坡[3]。然而根据对我国近100例发生较大变形或滑坡的斜坡进行的统计，其中反倾向边坡占33%[4]，出现的频度和危害较大[5]。块状-弯曲倾倒是反倾层状岩质边坡的主要倾倒形式之一，主要特性表现为坡顶后的弯曲倾倒以及坡顶前的弯曲倾倒的结合[6]。这种破坏模式在软硬相间的层状岩体中非常普遍，多发于砂岩泥板岩互层、燧石岩页岩互层、薄层状石灰岩中[7]。考虑到软岩、硬岩有不同的地质及工程特性，离心机模型试验和数值模拟成为研究此类坡体变形的主要研究手段。

关于离心机试验，Adhikary等[8]开展了7组不同结构的均质反倾边坡离心试验，得出了岩层倾角、岩石抗拉强度等对破裂面形状的影响；郑达等[9]在开挖条件下对均质层状反倾边坡进行离心试验，得出临空条件是深层倾倒变形破坏的关键致灾因子。上述成果均建立在相邻岩层岩性相同的假定条件之下，未考虑软硬相间互层状结构。

在数值模拟方面，黄波林等[10]使用UDEC软件对廖家坪危岩体进行了模拟，结果表明坡体倾倒的方向受软弱层控制；黄润秋等[11]以皖南某高速公路反倾层状边坡为例，使用离散元法分析了其倾倒破坏机理并揭示了边坡的变形分区现象。以上研究总结了坡体的破坏特征，但对于小变形状态下的细节兼顾不足；对于工程实例模型往往需要大量的标定来校核其微观参数，在同一参数不同边界条件下的继承性也略显不足。

白洁等[12]以苗尾水电站为背景构建反倾层状岩质边坡有限元模型，探讨了此类坡体的变形特点及稳定性影响因素；姚文敏等[13]采用FLAC3D强度折减法，在岩层倾角、岩层与边坡走向夹角变化时，探讨了三维软硬互层边坡的稳定性情况。以上研究由于方法限制无法进一步得到坡体的最终破坏情况。王宵等[14]以某梯级电站厂后边坡为例，分别构建二维离散元和三维有限元模型，对“似层状”岩质边坡具体演化过程进行了详尽的描述，但两模型的衔接方面并不一致，与实际条件不符。

F-DEM（有限元-离散元）方法[15]可以同时考虑岩体结构面和岩块，融合了有限元和离散元各自的优势，可以考虑岩层没有完全脱离母岩的实际情况。借助此方法，陈小婷等[16]以巫山县长江左岸危岩体为例，证明了将连续-离散模型用于边坡破坏研究具有可行性；刘郴玲等[17]以红石边坡工程为实例，建立FDEM模型并获得边坡失稳破坏的渐进破坏全过程。

本文首先对岩层倾角为优势倾倒角[18]（60°）的软硬互层反倾边坡进行试验研究，提出一种可以考虑拉剪和压剪的节理单元强度准则，并与软硬互层离心机模型试验相结合，分析了软硬互层岩质边坡中岩层倾角、软/硬岩层厚比对模型破坏机理的影响。

1 软硬互层倾倒变形体离心试验

1.1 试验设计

1.1.1 相似比设计

试验设计各模型/原型的比值如下：重力相似比（Ca）为90∶1、几何相似比（Cl）为1∶90、弹性模量相似比（CE）为1∶1、材料密度相似比（Cρ）为1∶1，黏聚力相似比（Cc）为1∶1，内摩擦角相似比（Cφ）为1∶1，应变相似比（Cε）为1∶1。最终预堆砌模型边坡的尺寸为70.27 cm（长）× 55 cm（宽）× 56.72 cm（高）。

1.1.2 相似材料设计

石膏水泥混合材料是目前模拟岩石行为最常用的相似材料，可用于控制材料的力学性质。为避免硬化过程中受环境湿度与温度的影响，加入硼砂水（m硼砂∶m水=1∶49）。经正交配比试验，再通过单轴压缩试验和三轴压缩试验，得到原型和相似材料的物理参数（表1），最终材料配比见表2。通过100 mm×50 mm试样单轴压缩试验，得出软岩相似材料抗压强度为0.3 MPa，硬岩为10 MPa。

表1 原型材料与模型材料的物理力学参数Table 1 Physico-mechanical parameters of the prototype andmodel materials

表2 软硬岩相似材料最终配比Table 2 Ratio of similar materials for hard and soft rocks

1.2 试验设备与模型加载

1.2.1 模型信息

本文软硬互层边坡离心物理模型试验在地质灾害防治与地质环境保护国家重点试验室（成都理工大学）TLJ－500型岩土离心机上完成。

边坡模型岩层倾角为60°（优势角[18]），由单块板（岩层）拼装而成，试块分别有60 cm×10 cm×1 cm和40 cm×10 cm×1 cm两种类型。如图1所示，涂黑岩层为硬岩层，未涂黑岩层为软岩层，岩层厚度均为1 cm。为避免边界效应和降低侧壁摩阻力的影响，在模型箱与模型接触的位置粘贴光滑的塑料薄膜。

图1 模型及监测点信息Fig.1 Model and monitoring point information

1.2.2 岩体节理及层间处理方法

岩层搭接方式为错缝搭接（图2），且相邻两块板使用软石膏砌筑加强黏结，这样搭接可以大限度避免垂直于坡面方向节理的影响[9]。

图2 模型搭接方式及节理情况Fig.2 Overlap mode and joint of the model

在层状坡体的制作过程中，相邻层理面宜界面均匀、层厚均匀，所以在浇筑拼装时须注意控制时间，每浇筑完成一层，养护1 h再浇筑下一层。层间黏结力宜较小，因此层间材料选用抗压强度0.1 MPa、抗拉强度10 kPa的软石膏。

1.2.3 加载方式

由离心机基本原理可知：

式中：am—离心加速度/（m·s−2）；

N—离心机转速；

g—重力加速度/（m·s−2）。

离心机模型的尺寸将由离心加速度放大至N倍，由此使较大尺寸的边坡变形破坏过程的模拟试验成为可能。

模型所受离心加速度首先从0g开始逐渐加载至30g，并在30g保持离心机匀速旋转5 min；然后加速度从30g逐渐加载至50g，并在50g保持匀速旋转5 min，最后模型所受加速度从50g逐渐加载至90g，并在90g保持匀速旋转5 min。若模型不破坏则继续增加离心加速度直到模型破坏。

1.2.4 监测点布置

本研究采用导轨支撑位移计量测坡体表面位移（图1）。在坡表设置两个监测点：监测点A和监测点B，并在层间埋置压力传感器监测层间法向应力。

2 软硬互层倾倒变形体数值模型

2.1 F-DEM模型介绍

F-DEM模型由实体单元和零厚度的节理单元两种单元构成（图3）。其中，零厚度的节理单元夹在实体单元中间。当节理单元受力过大时，程序在当前增量步完成后将节理单元删除，并在删除位置生成对应的裂缝。节理单元删除后的两个实体单元间依旧保持接触关系，可以碰撞和摩擦，法向接触为硬接触，切向为库伦摩擦，其中硬接触如图4所示。

图3 F-DEM模型示意图Fig.3 Schematic diagram of the F-DEM model

图4 法向硬接触Fig.4 Normal hard contact

对于三维单元，有

式中：E—弹性模量/MPa；

μ—泊松比；

σx、σy、σz—x、y、z方向应力；

εx、εy、εz—x、y、z方向应变；

γzx、γxy、γyz—x、y、z方向切应变。

本文所使用的节理单元假设剪应力为0（在局部可以将沿着节理方向的两个主应力看作剪应力）[19]，所以式（1）可简化为：

由于剪应力都为0，所以在节理单元中只存在3个主应力：σ1、σ2和σ3。在局部坐标系（图5）中：

图5 局部坐标系Fig.5 Local coordinate system

2.2 Hoek-Brown与Mohr-Coulom联合强度准则及VUMAT子程序实现

Hoek-Brown准则和Mohr-Coulomb准则是常用的用来描述岩石抗剪强度随法向应力变化规律的强度准则。Hoek-Brown准则是根据三轴压缩试验数据得到的经验公式：

式中：σ1、σ3—最大和最小主应力/MPa；

σ0—岩石的单轴抗压强度/MPa；

m—与岩性相关的常数；

S—岩石完整性的参数，完整岩石时S= 1。

另外，岩石的单轴抗拉强度可用式（5）预测：

Mohr-Coulomb准则用式（6）表示：

Hoek-Brown准则比Mohr-Coulomb准则能够更好地描述岩石在拉剪应力状态下的强度特征，而Mohr-Coulomb准则比Hoek-Brown准则能够更好地描述岩石在压剪应力状态下的强度特征[20]。采用分段函数的方式来对砂岩拉剪强度和压剪强度进行预测：

其中式（4）中的第二项又可写为：

因此将Hoek-Brown准则与Mohr-Coulomb准则联合使用，既能较为准确地描述岩石在拉剪应力状态下的强度特征，又能较为合理地揭示岩石在压剪应力状态下的强度特征，可以更好地描述（模拟）岩石的断裂行为。VUMAT子程序的计算目的是联合两种强度准则以及生成新的裂缝。每一次迭代过程都将遍历模型中各个单元。若判断为受压状态，则进一步判断其当前强度是否超过其Mohr-Coulomb准则下抗剪强度，若大于抗剪强度则将损伤因子D更新为1，清零其刚度矩阵。此时单元失效。若判断为受压状态，则进一步判断其当前强度是否超过其Hoek-Brown强度准则下的抗剪强度，若大于抗剪强度则将损伤因子D更新为1，清零其刚度矩阵。此时单元失效，裂缝产生。此时实体单元之间的接触关系依然存在，与母体分离后的岩块依然可以相互摩擦、碰撞。其中σn>0时模型处于压剪状态，σn≤0时模型处于拉剪状态。

2.3 F-DEM模型构建

2.3.1 F-DEM中的时间及零厚度单元插入

使用ABAQUS中Dynamic、Explicit分析步进行计算，此算法中的时间即为运动方程中的时间，是真实时间。编写Python程序对有限元模型进行前处理，步骤如下：

（1）将原始模型中的节点分别赋予两个不同的序号，其中一个节点为虚拟节点；

（2）将新产生的节点编号以及节理单元节点的编号添加到模型文件中；

（3）添加节理单元的数据信息，即在源文件的关键字中声明单元类型的更改；

（4）为这些新加入的单元设置材料和截面属性并再次导入ABAQUS中进行边界条件等的设定。

2.3.2 Hoek-Brown—Mohr-Coulomb联合强度准则和本构模型参数

图6为软硬互层反倾层状边坡数值模型图，模型按单元不同分为两部分：实体单元部分和零厚度的节理单元部分。

图6 软硬互层反倾层状边坡数值模拟Fig.6 Numerical model of anti-dumping rock slope interbeded by hard and soft layers

根据坡体岩性和结构面，将节理单元分为层间节理单元、硬岩节理单元及软岩节理单元。因VUMAT中的参数都为实际物理量，所以可以依照坡体离心机试验模型中的材料参数设置数值模型中的参数，且模型建立两者一致（其中Hoek-Brown模型中的参数经过换算得到[21]）。对于实体单元，则直接使用单轴压缩试验所测得的材料参数。对于层间节理单元的参数通过与试验现象的对比校核来确定。通过以上材料参数的确定方法来确保物理模型与数值模型两者之间材料的相似。数值模型各部分材料具体材料参数见表3、表4。

表3 节理单元参数Table 3 Joint element parameters

表4 实体单元参数Table 4 Solid element parameters

2.3.3 计算工况及条件设定

将10 120个零厚度节理单元嵌入到6 828个实体单元之间。为确保物理模型与数值模型两者之间边界条件一致，离心机模型（图1）按1∶1比例构建相应F-DEM模型，模型尺寸为700 mm×570 mm，约束所有节点厚度方向的位移。为保证数值计算在加载条件上与离心机试验一致，数值模型按1.2.3节中所述加载条件进行增重加载。

3 离心试验结果与模拟结果对比

3.1 监测点位移对比及法向应力监测

监测位移曲线及数值计算得出的位移曲线对比情况见图7（监测点布置见图1），二者变化趋势高度接近，初步说明了所建立数值模型的正确性。

图7 监测点位移对比Fig.7 Monitoring point displacement

利用压力传感器对模型的层间法向应力进行监测，得到监测结果（图8）。可知，层间法向压力监测点b1、b2、b4的压力较早上升且伴有波动，加载重力进一步增大至50g时，边坡b1、b2监测点的压力值进一步增加。当加载接近90g时，随着深层弯折带的贯通，压力监测点b1、b2的值到达最大。

图8 层间法向压力监测曲线Fig.8 Monitoring of normal pressure between layers

3.2 破坏过程对比及数值模型正确性验证

图9（a）为离心加速度10g时坡体的变形情况。根据图8的位移曲线，加载至20g之前，边坡处于稳定阶段，坡内无可见变形；此时监测点b1、b2所测得的应力也未到达第一次峰值。超过20g后，位移量出现逐渐增大，持续至30g，后缘出现拉张裂缝，坡顶及坡表出现轻微的相互错动，此时监测点b1、b2、b4也出现应力波动。加载至30g时，位移曲线出现小幅跃升（图7），此时坡体已处于破坏临界状态。图9（b）为此时变形情况，坡脚局部岩体折断、裂缝明显，且坡顶拉张裂缝数量迅速增多。50g时，位移急剧增大，边坡变形如图9（c）。坡体内一级破裂面已接近贯通。与此同时b1、b2监测点的应力迅速增大，因为此时破裂面以上岩体发生倾倒变形，坡脚岩体被挤出，形成脱空，加剧了上部岩体的变形。此时被挤出岩体的末端已开始形成更陡的二级破裂面。90g时，监测点A、B的位移以及监测点b1、b2的应力达到峰值，并暂时维持稳定，此时边坡变形见图9（d），第二破裂面已向上贯通至坡顶，并直接导致上部岩体倾倒变形，产生三级破裂面。此时几乎所有层间法向压力计的测量值都有所增加，说明层间空隙被层间法向压力挤实，裂缝的发展接近稳定状态，变形破坏达到最终形态。

图9 数值模拟与试验对比Fig.9 Comparison between numerical simulation and experiment

通过以上分析可知，离心机试验模型的变形破坏情况与F-DEM数值计算模型基本一致，验证了数值模型及所提出强度准则的正确性，且层间应力与于此类边坡变形破坏的关联性较大。为能够灵活模拟不同层间强度因素的影响，建议将层间节理单元的本构修改为“牵引-分离”准则。

4 不同几何结构特征对边坡倾倒特征的影响

4.1 数值试验设计

岩层倾角以及软/硬岩层厚度比是反倾软硬互层边坡的主要几何因素[11]，坡顶位移以及破裂面形态可以表征倾倒破坏特征。固定软/硬岩层厚比为1∶1（10 mm∶10 mm），选取40°、60°、80°岩层倾角，确定方案1用以对比分析不用岩层倾角对边坡倾倒破坏特征的影响；固定岩层倾角为60°选取软/硬厚度比1.5∶1（15 mm∶10 mm）、2∶1（20 mm∶10 mm）、1∶1（10 mm∶10 mm）、1∶1.5（10 mm∶15 mm）以及1∶2（10 mm∶20 mm）确定方案2，用以对比岩层软/硬岩层厚度比不同对边坡倾倒破坏特征的影响。

为保证数值试验方案的关联性与合理性，所建不同的倾角和不同软/硬岩层厚度比的模型均按照2.3.3小节中的边界条件进行设置。

4.2 岩层倾角对边坡倾倒破坏的影响

4.2.1 岩层倾角对边坡最终破坏形态的影响

岩层倾角40°时，边坡达到最终破坏时的离心加速度为100g，坡体后边局部小范围张开，一级破裂面与水平面的夹角（锐角）近似等于60°，且坡脚处有较小范围的岩块剪出，如图10（a）所示。岩层倾角60°时，坡体共产生两条贯通的破裂面，坡体沿二级破裂面滑出，如图10（b）所示。相对于岩层倾角40°时边坡的最终破坏形态，岩层倾角60°时坡体后缘的张开范围以及坡脚剪出范围更大，且一级破裂面的分布位置更深。此时边坡达到最终破坏时的离心加速度为90g，相对于岩层倾角40°、60°时，岩层倾角80°时，边坡达到最终破坏时的离心加速度为105g，边坡更倾向于整体破坏，后缘张开范围较大，如图10（c）所示。这说明随着岩层倾角的增大，破裂面由一条逐渐向多条发展；且岩层倾角越大，边坡破坏的整体性越强。

图10 不同倾角下边坡最终破坏形态Fig.10 Final failure pattern of slope under different inclination angles

由以上分析可知，岩层倾角对第一破裂面的影响较为显著（图11），岩层倾角40°时，边坡一级破裂面较为陡峭且呈现出连续的具有台阶状的滑移面。随着岩层倾角的进一步增大，破裂面向坡体更深处发展，“台阶状”形式减弱。随着岩层倾角的增大，边坡一级破裂面逐渐向坡体深处发展，当岩层倾角80°时，边坡一级破裂面发育位置较深，且后缘拉裂缝较大，边坡发生整体倾倒，说明此时边坡发生以“弯曲（倾倒）-拉裂”为主的破坏[3]。

图11 不同倾角下边坡一级破坏面Fig.11 First-order fracture surface with different dip angles

4.2.2 岩层倾角对坡体位移特征的影响

图12为监测点A所记录的位移-加速度曲线。在滑坡启动阶段，重力加速度较小，不同倾角对监测点A的竖向位移影响不大。随着重力加速度的增大，倾角越大，坡体初始滑移所需重力加速度越大。岩层倾角越大岩层垂向的重力分量越小，所以倾角越大初始滑移需要的重力加速度越大，进一步增重后进入最终阶段形成滑坡。岩层倾角越大最终产生位移越大，坡体破坏程度越剧烈。在启动阶段达到初始滑移所需的重力加速度越大。重力做功下，坡体积累的总的弹性应变能越大，突然释放时所产生的动能也就越大。所以，岩层倾角越大最终阶段产生的位移越大。

图12 不同岩层倾角下竖直位移Fig.12 Vertical displacement with different dip angles

4.3 软/硬岩层厚度比对边坡破坏特征的影响

4.3.1 软/硬岩层厚度比对边坡最终破坏特征的影响

为探究软/硬岩层厚度比对边坡破坏特征的影响，将倾角60°，离心加速度90g时，不同软/硬岩层厚度比的边坡破裂面绘制于图13。如图13（a）—（d）所示，坡体后缘存在两处折断区域。随着软/硬岩层厚度比由1∶1.5向2∶1转变，两个折断区之间距离越来越近，最终合为一个折断区。随着硬/软岩层厚之比增大，坡体后缘折断区深度逐渐增加。坡体在增重过程中的附加应力主要由硬岩承担。硬岩岩层所占比例越大，坡体增重过程中所积累弹性应变能越大，坡体破坏释放的能量就越大，坡体滑动的整体性就更强。坡体滑动的整体性也表现在破裂面的分布深度上，对比图13（a）—（e）可知，二级破裂面的分布深度随着软/硬岩层厚比的提高逐渐减小。

图13 不同软/硬岩厚度比边坡最终破坏形态Fig.13 Final failure pattern of slope with different ratio of soft/hard thickness

硬岩岩层厚度越大，能承受的弯矩越大，则二级破裂面分布越深。破裂面的形态也发生了巨大变化：随着软/硬岩层厚度比的增加，二级破裂面的形态逐渐由粗糙的“锯齿状”向平滑的“圆弧状”转变。这是因为软岩层厚所占比例越大，越接近土质滑坡的大弧度破裂面。坡体中的一级破裂面也有类似特征。

4.3.2 软/硬岩层厚度比对边坡位移特征的影响

如图14所示，随着加速度的增加，坡顶的竖向位移呈阶梯式增加。软/硬岩层厚度比为2∶1时，坡顶竖向位移最大（约250 mm）。随着硬岩层厚所占比例的提高，坡顶竖向位移逐渐减小；软/硬岩层厚度比为1∶2时，坡顶竖向位移最小（约170 mm），这是因为硬/软岩层厚度比值越大，坡体的整体弹性模量越大。

图14 不同软/硬岩层厚度比下坡顶竖直位移Fig.14 Vertical displacement of downhill top with different ratios of soft/hard thicknesses

5 结论

（1）经离心试验验证，在F-DEM模型节理单元中使用Hoek-Brow与Mohr-Coulomb联合强度准则VUMAT子程序可以较好地模拟软硬互层反倾层状岩质边坡破裂面的大致开裂情况。

（2）软/硬岩层厚度比为1∶1，岩层倾角60°的坡体破裂面大致有三个。其破坏过程为：坡脚部位开始出现弯曲变形，边坡后缘出现拉张裂缝，边坡整体向临空面弯曲倾倒。从坡体内部至坡表依次形成三个破裂面。

（3）随着岩层倾角的增大，边坡一级破裂面逐渐向坡体深处发展，且坡体由“弯曲-拉裂”为主的破坏模式向“滑移-压致拉裂”为主的破坏模式转变。岩层倾角越大最终产生的位移越大，坡体最终的受破坏程度也就越大。

（4）软岩所占比例越大坡顶位移越大；硬岩所占比例越大坡顶位移越小。且不同比例的软/硬岩层厚度对坡体滑动的整体性以及破裂面形态有一定影响：硬岩所占比例越大，坡体滑动的整体性越强；软岩所占比例越大，坡体的破裂面越接近土质滑坡的“圆弧状”。