基于弹塑性有限元的深挖路堑碎裂岩质边坡稳定性分析

秦 浩，陈千寻，张 华

（1.西华大学流体及动力机械教育部重点实验室，四川 成都 610039；2.西华大学能源与动力工程学院，四川 成都 610039；3.成都大学建筑与土木工程学院，四川 成都 610106）

随着我国经济的飞速发展和西部大开发，西部地区高等级公路相继开工建设。西部地区地形地质条件复杂，前期的工程勘察和试验工作有限，对工程对象的复杂地质条件始终难以把握[1]。边坡工程需要因地制宜地进行设计，工程地质条件和施工设计的不同促成了边坡的特殊性，对于具有复杂地质条件的边坡，利用以往经验对其进行常规的稳定性分析不能得出合理充分的结论。因此，如何对岩质边坡进行有效合理的稳定性分析已成为西部地区工程建设中的重要问题。

边坡工程的稳定性受控于其岩土体的基本特性和人为改造程度。岩土体具有多变性和不均质性，随着施工的进行会发生一系列动态变化影响边坡稳定性，在对边坡进行设计时需要依据已有信息预测其发展趋势。信息化施工利用施工中逐步揭示的工程地址信息，结合数学、力学分析方法对施工过程应力、应变进行计算分析，并依据分析结果对工程进行及时设计调整以达到最优效果[2 − 4]。在信息化设计的基础上，将实际工程中的监测和理论分析结合起来，对工程岩体的变形进行预测性分析，为施工方案的优化设计及边坡稳定性分析提供依据[5]。

传统的边坡稳定分析主要使用极限平衡法，但极限平衡法会对问题进行简化与假定，与实际情况之间有不同程度的差异[6]。在常用的数值计算方法中，离散元法能考虑块体旋转、脱离等情况[7]，在处理散体结构的“散”“动”特性上有较好的效果，但计算精度并不比有限元法高很多。而有限元法能从区域整体上掌握岩体受力、变形规律，适用于各种复杂形状，在大规模工程稳定性分析上具有离散元法不可比拟的计算效率优势[8]。在有限元法的研究上，Griffithe 等[9]通过有限元法研究了边坡的变形特征、破坏特征和稳定性。郑颖人等[10]对有限元强度折减法在边坡分析中的破坏判据、计算精度和影响因素进行了研究，并用此法对边坡进行了稳定性分析。杨建平等[11]利用有限元方法模拟工程裂隙岩体加载破坏过程，对计算岩体强度参数取值进行了研究。

虽然数值模拟技术在边坡稳定性分析中给出了较高的精度，得到广泛应用，但是在计算过程中需要解决的包括地质力学模型、力学参数等关键问题上还无法给出准确的答案[12]。BP 神经网络优良的数据拟合能力和建模能力在非线性数值计算中表现出较好的适用性，使用BP 神经网络反演岩体力学参数能获得令人满意的结果[13]。

本文针对碎裂岩质边坡，采用将有限元法和BP 神经网络相耦合的方法，利用施工现场监测的位移值反算有限元模型参数，获得符合实际情况的最优岩体参数，并利用算得的模型参数对边坡稳定性和开挖加固效应进行预测。本研究为信息化施工设计提供了一种新方法。

1 工程概述

1.1 地质概况及防护设计

该路堑边坡工程地表多被第四系崩坡积层掩盖，斜坡平均坡度约45°，顺坡自上而下分为砂质泥岩强风化区、泥质砂岩中风化区、厚度约1 m 的凝灰岩中风化区与软弱夹层，中间包含一个平均宽度约10 m 的断层带。受断层影响，该边坡中岩体极为碎裂，风化严重，岩层裂隙发育且倾角较陡。边坡岩层整体产状较为稳定，节理面以陡倾为主并与路线斜交，岩层面陡倾且垂直于坡面（见图1）。

图1 施工边坡地质剖面图

1.2 施工监测设计

施工监测采用地表监测和地下监测相结合的立体监测系统，以位移监测为主。

在坡顶和平台设置地表变形监测设施，监测地表变形。综合地形地质条件，采用测斜的方式监测边坡工程内部的水平变形，在边坡第一、第三、第五级平台设深部位移监测孔，3#监测孔设于第三级平台上，孔深29 m，5#监测孔设于第五级平台上，孔深27 m。

2 有限元模拟计算与参数分析

边坡具有非线性、非均质特性，使用有限元法进行模拟计算边坡位移场和应力场变化，能够较接近实际地反映边坡的破坏机制。依据工程实际情况建立有限元求解模型，将原本的连续域离散化为有限个单元，对单元结点的状态变量进行求解，并选取合适的准则对求得的解进行分析，判断边坡状态及发展趋势。

2.1 计算原理

岩土材料同时具有弹性和塑性性质，边坡稳定性分析是对边坡变形甚至破坏的预测性计算分析，而在其变形过程中的弹性和塑性两个阶段是密切联系的，应用弹塑性模型为有限元模型的本构模型。弹塑性模型计算分析是使对象整体处于弹性阶段，允许局部区域进入塑性状态，由于塑性变形具有不可恢复性，因此采用塑性力学破坏理论描述其模型塑性变形及破坏特征[10]。相应地，采用岩土材料常用的德鲁克−普拉格准则（D-P 准则）判断岩体从弹性到塑性的状态变化能获得可靠性较高的计算结果[14]。岩石的力学行为和参数在加荷或卸载时是不同的，可以将开挖看作是对岩石施加的一种荷载，因此整个开挖过程对每一个特定部位等同于一个反复加卸载过程[15]。对开挖边坡的稳定性进行分析，需要对每一步开挖后边坡应力场变化进行模拟，采用施加等效释放荷载分析计算受开挖影响的回弹变形量，等效释放荷载计算的正确性直接影响稳定性分析的结果。

采用增量理论描述边坡在施工过程中变形增量与应力、应变增量的关系。岩体进入屈服后，边坡模型的全应变增量d{ε}可表示为弹性应变增量d{εe}和塑性应变增量d{εp}之和：

其中：[Dep]为与边坡单元材料相关的弹塑性矩阵，由相应的弹性矩阵 [D]和塑性矩阵[Dp]作差求得，[Dep]=[D]−[Dp]；{σ}为边坡单元内任一点的应力分量。弹性变形服从Hook 定律：

式中：G为外部荷载进入岩体弹塑性变形区域时岩体的塑性势函数，G=G({σ},Ha)，Ha为岩体进入屈服处于弹塑{性}变形状态产生应变硬化的硬化参量，λ为塑性势函数G对应力{σ}偏导得到的应变量与塑性应变增量的比例系数，λ≥0。

弹塑性应力应变关系的一般表达式为

相比Tresca 准则和Mises 准则，D-P 准则同时考虑了静水压力与中间主应力对岩土体屈服强度的影响，更加适用于岩土材料。复杂应力下采用Mohr-Coulomb 外角点外切圆准则对应的D-P 准则判断岩体状态可获得较大塑性区。

D-P 准则通过屈服函数f表示为：

式中：I1为岩土材料应力张量的第一不变量；J2为岩土材料应变偏量的第二不变量；α，K为关于岩土体强度参数的常数。它们采用下式确定：

其中：φ为岩土材料摩擦角；C为黏聚力。

当岩体从屈服进入塑性状态，此时岩体应力满足屈服函数值f=0，取与屈服函数表达式相同的塑性势函数，塑性矩阵为：

其中，A为应变硬化参数。

在施工过程中开挖会使部分原始边坡内部单元面暴露在空气中，形成一系列边界点，在施工过程中各点应力在加载或卸载间反复变化，开挖后边界点的应力被解除使得岩体原始应力场发生改变。挖去单元e 使得原本处于岩体内部的单元面成为自由表面，在计算中要将单元e 在原始应力场中对相应内部单元面的作用力等效加载到相应自由面上。

自然状态下岩体的应力场处于平衡状态，当挖去一个单元e 时，其相邻结点i受到的原始应力条件被改变，对i点施加单元e 产生的结点力Fi和等效荷载Xi形成新的平衡应力场。

其中，[B]为单元应变矩阵；{σ0}为开挖前单元e 的应力分量；Ωe为挖去单元e 的体积。

如果单元e 被开挖掉，对i结点产生应力释放，其释放力Ri为：

假设被开挖掉的单元为k个，则总的释放荷载引起的等效结点荷载{R}为：

式中：[Ce]为单元结点力加载到单元体上的系数矩阵；γ为岩体沿各坐标方向的容重分量所构成的列阵，通常情况下，地下岩体所受荷载仅是它的自重。

2.2 有限元计算模型

边坡工程的计算模型要反映边坡各种结构的空间几何关系和拓扑结构关系，将各种结构面和边界面准确划分出来，满足工程复杂性和多样性，缩小其不确定性。基于对施工边坡的工程地质结构的解析，根据其地层岩性和地质构造的几何形态和空间关系将模型分为7 个部分，如图2 所示，其中1和6 为泥质砂岩强风化区，2 到5 分别为泥质砂岩中风化区、软弱夹层、凝灰岩中风化区、断层带。

图2 边坡材料分区示意图

计算范围取垂直路线方向宽(X向)115 m，高度(Y向)取90.5 m，顺公路方向长(Z向)82 m。右边界取到第五级坡面顶部向山体内延伸50 m 处，左边界由地形图确定。

网格划分的好坏直接关系到有限元分析计算的精度和速度。根据边坡的工程地质特性，从计算的精度、时间和建模工作量等方面综合考虑，依据各分块区的特点采用20 节点六面体实体单元及其退化形式的混合网格进行有限元模拟。对于边坡底部结构性质稳定，受开挖影响相对较小，采用精度较高的六面体网格，而表层风化严重、岩体破碎，采用四面体网格。

边坡模型的边界条件设为底部固定，在模型邻公路面和背面施加水平向双向约束，在垂直公路两侧施加法向约束，坡体表面设为自由面。因为主要在表层碎裂岩体开挖，只涉及边坡表层岩体卸荷，所以将有限元开挖模拟计算的初始应力场设为自重应力场。未开挖状态下的三维边坡模型网格划分如图3 所示。

图3 边坡网格划分结果图

2.3 参数分析

岩体的岩石力学参数通过对现场钻孔安装测斜管取回的岩样进行岩石力学试验初步确定。采用烘干试样进行常规三轴试验获得不同围压下岩石的抗压缩强度，结合单轴抗拉、抗压强度试验计算得出内摩擦角φ和黏聚力C。各分块岩体力学参数初步取值结果见表1。岩样的获取是随机性的，只能代表相应局部范围的岩体物理性质，并且从取样到试验这一时段有大量的外在因素会对岩样物理性质造成影响；因此，通过岩石力学试验获得的岩体参数直接用作有限元模型的计算参数会导致计算分析结果失真。但是通过岩石力学试验可以大致确定岩体参数的范围，其初步取值结果能作为基础参数数据，为进一步计算参数的准确取值分析计算研究提供支持。

表1 岩体力学参数初值

3 参数反分析

由于岩石试样采集的局限性和岩石力学试验的有限性，各区材料参数的选取不可避免地与实际情况有所偏离。岩体弹塑性模型力学参数包括密度ρ、弹性模量E、黏聚力C等对模型位移计算结果的影响是不同的，并不是所有的参数都需要参与参数反演，在反算前获得各参数的位移影响度大小能确定其可反演性，辅助计算参数反分析的进行[16]。本文利用BP 人工神经网络建立岩体参数与位移的非线性关系，通过参数敏感性分析来扩大参数样本数量提供给神经网络，对比神经网络的位移计算值和实际监测位移值，不断缩小两者差距获得逼近于实际的模型参数。

3.1 敏感性分析

计算方案为模拟从三级平台开挖到路面这一施工周期的岩体位移变化。敏感性分析采取对每个分块区的每个岩体参数依次进行折减，在对任一个参数进行折减时保持其它参数不变，以折减计算所得的位移值与初始计算所得的位移值的差值大小判定模型对各参数的敏感性[17]。折减方案如表2所示，在岩石力学试验获得的岩体力学参数基础上，方案实施顺序按表1 块体序号从左到右、从上到下依次对表内参数进行折减，每个方案只折减一个岩体块区的一个参数。将每个方案的参数代入弹塑性有限元模型进行计算[18]，各方案顺坡向位移变形计算结果如图4 所示。

表2 参数折减方案

图4 各方案顺坡向变形计算结果

由图4 可知，与未折减方案相比，方案2、5、6、10、14、19、22、26 表现为顺坡向变形减小(朝向坡内)，方案1、7、11、15、18、23、27 表现为顺坡向变形增大(指向坡外)。其余折减方案计算结果显示顺坡向变形很小。进一步分析可得，当1、5 区的密度ρ折减，或是2、3、4、6 区的弹性模量E折减时，坡体呈现顺坡向变形增大的趋势，不利于边坡稳定；当2、3、4、6 区的密度ρ折减，或是1、5 区的弹性模量E折减时，坡体呈现顺坡向变形减小的趋势。还可以看出，内摩擦角φ和黏聚力C的所有折减方案的计算结果都表现出顺坡向变形微小，参数敏感性弱[19]。

另一方面，分区1 为表层碎裂岩体也是边坡主要开挖区域，而开挖必然会使得1 区的弹性模量E和黏聚力C发生改变。分区3 为软弱夹层，处于中间层(分区2)与基岩(分区4)之间，2、3、4 区对各参数的折减反映出相同量级的发展趋势，当弹性模量E折减时，3 个块区都表现出顺坡向位移增大的趋势，并且表现出越靠近坡面区顺向位移变化越大的现象。分区5 为断层带，岩体极为破碎，但是对其弹性模量E进行折减并不会使得边坡向不稳定方向发展。从对内摩擦角φ和黏聚力C的计算结果看出，这两个参数对边坡的变形影响较弱，反应出施工边坡抗剪强度高。

由参数敏感性分析结果得出，要使得施工边坡稳定，可通过改变其密度ρ和弹性模量E来实现。而工程对象属于岩质边坡，通过改变岩体密度来实现边坡稳定对于实际工程来说可行性低，而改变弹性模量E可通过预桩锚固实现。

3.2 模型参数反分析

建立反映边坡位移变形量与计算参数非线性关系的BP 人工神经网络模型，将实测位移值和样本集里的位移值代入BP 人工神经网络模型进行参数反算，将样本参数代入有限元模型计算位移值，通过多次反向和正向计算分析，不断逼近实际监测位移值，获得最优有限元模型计算参数，为后续边坡稳定性分析提供依据[20]。

将经过敏感性分析扩大的参数样本作为BP人工神经网络的训练样本，当输出误差满足要求则训练结束。将实测位移值代入已训练好的BP 神经网络模型对各岩体参数进行预估，再将预估得到的密度ρ、泊松比μ、弹性模量E、内摩擦角φ、黏聚力C代入有限元模型中进行开挖模拟计算得到相应计算位移值，通过多次BP 神经网络和有限元法耦合计算，在参数反演和正向分析中不断缩小与实测位移值的差距，最终获得最优的模型参数。经由反算得出的各岩体参数见表3。

表3 经过反算得出的参数表

由表1 和表3 可知，弹性模量E、内摩擦角φ和黏聚力C在进行反算前后在每个块区变化明显，密度ρ没有变化。密度ρ经过反算后，每个分块区的值和通过岩石力学试验确定的初值相同，其中，5 区的密度ρ最小，为2.5 g/cm3；7 区最大，为2.78 g/cm3，4 区次之，为2.6 g/cm3。弹性模量E反算前和反算后的数值有着相当大的差距且整体呈下降趋势，在其数值大小上，反算前后的最小值都出现在5 区，最大值都出现在7 区，次大值在4 区；在其变化量上，1 区的变化量最大，反算后的弹性模量是反算前弹性模量的54 倍，2～7 区反算后的弹性模量是反算前弹性模量的10 倍。内摩擦角φ反算后比之反算前降低，3 区下降5°，其余各区均下降了10°。黏聚力C反算后比之反算前，总体呈下降趋势，在1 区下降了94%，2 区下降90%，3 区下降67%，4、5、6、7 区分别下降了59%、50%、83%、87.5%。

对比表1 和表3 各岩体参数反算前后的数值变化可得：密度ρ在反算后整体表现稳定，在各区数值大小的不同反映出了不同岩质块区的密实度的不同。弹性模量E与反算前相比显著减小，在模拟开挖卸荷过程的有限元参数中表现出相当大的变化性，1 区为边坡主要开挖的表层碎裂岩体区，其弹性模量E缩小最多，3、5 区分别为软弱夹层、断层带，但其受工程开挖的影响与其余岩质较好的块区相同。在内摩擦角φ的变化上，模拟开挖过程反算后，整体呈下降趋势，但2、4 区之间的3 区软弱夹层相比其他各区减少量小，5 区断层带内摩擦角φ与其紧邻的2、4、6、7 区减少量相同。在黏聚力C的变化上，模拟从一级开挖三级平台过程使得岩体中的含水量随着与空气接触量的增加而下降，表现为开挖区1 区C的下降最多，紧邻开挖区的2 区次之，软弱夹层3 区和断层带5 区下降最少。反映出了软弱夹层和断层带的强度稳定性良好，对边坡的变形破坏影响小，而表层开挖区1 区及其紧邻区2 区受开挖直接影响较大。

图5 为利用反算得到的岩体参数模拟开挖过程有限元位移计算结果和实际位移测量结果曲线图，图中结果在走势、形状、突变位移大小上都表现出显著的一致性，表明分析所得的岩体力学参数能够表现实际岩体的物理性质，具有准确性，用于边坡稳定性分析预测能够很好地贴合实际。

图5 计算的结果和实际测量结果图

另一方面，也反映出利用岩样进行岩石力学试验获取参数的局限性。各种工程措施下边坡的物理性质会发生必然的改变，岩样本身从岩体内取出到进行岩石力学试验的过程中也包含任意性和变化性，致使岩体参数具有很强的不确定性。耦合BP 人工神经网络和有限元模型对岩体参数进行反算，其本身就包含一个模拟边坡的动态变化过程和地质条件揭示过程，能反映出由于开挖造成的岩体含水量减少致使抗剪强度减小及开挖卸荷作用下岩体的回弹变形趋势。对比反算前的岩体参数，模拟开挖过程反算结果显示，岩体密度ρ不变，弹性模量E、内摩擦角φ和黏聚力C整体下降，结合敏感性分析结果可知，开挖会导致边坡有滑塌的趋势。模拟结果验证了复杂工程地质条件岩体力学参数智能反演的可行性[5]。

3.3 施工过程动态模拟

将分析所得的岩体力学参数应用于加固效应分析及边坡稳定预测。图6 为直接开挖顺坡向位移计算结果，开挖坡顶位移是负向的，但在坡脚处有一个正向的位移，整个边坡有滑塌的趋势。图7是在坡脚进行预加固方案的计算位移等值线图，结果显示进行预应力桩加固对边坡坡体产生的扰动和破坏较小，坡顶和坡脚的位移变形都朝向坡内，会引起边坡滑塌的坡脚顺坡向位移变形转化为对边坡稳定性影响相对较小的朝坡内的变形。对整个开挖坡面都采用锚索框架梁加固的计算结果见图8，根据设计施工方案，加锚力按单根锚索锚固力932.4 kN，方向与水平面呈20°夹角的加锚方式设计，此种方案同样能抑制边坡滑塌。

图6 直接开挖方案顺坡向位移等值线(mm)

图7 桩预加固方案开挖前顺坡向位移等值线(mm)

图8 锚索框架梁方案开挖后顺坡向位移等值线(mm)

从位移场的计算结果可看出，直接开挖方案，在边坡无防护状态下直接开挖对边坡本身稳定的应力场破坏严重，坡脚处产生顺坡向位移有致使边坡有滑塌的趋势，不利于边坡的稳定。对坡脚进行预加固方案的位移变形均指向坡内，表明坡脚预加固可有效防止坡脚应力集中现象和边坡后缘拉应力区的发生，抑制边坡滑塌。坡脚预加固方案的位移场计算结果还反映出边坡在原始自然状态下整体相对稳定，开挖卸荷作用使边坡原始应力场发生改变从而影响其稳定性。采用开挖坡面全锚索，使整个边坡呈现预应力受力状态，弹性模量E增大，坡体位移变形朝向坡内。对比坡脚预加固方案和开挖坡面全锚索框架梁方案，两种方案都能有效抑制滑坡，达到边坡稳定性要求，从经济上考虑预桩锚固方案相对较好。

4 结论

利用弹塑性有限元模型，结合信息化施工获取由开挖逐步揭示的地质信息，利用有限元模型模拟边坡开挖过程位移变形计算，对岩石力学试验获得的岩体力学参数初值进行敏感性分析，并进一步进行参数反分析，当位移计算值和实测值接近程度最大时，即获得最优模型参数，应用于边坡稳定性预测分析和加固效果预测，得到了合理的计算结果，并得出以下结论。

1）采用弹塑性有限元对边坡进行稳定性预测和加固效应分析是有效可行的，可以有效地计算得出岩体位移场、应力场的变化特征，进而得出有破坏趋势的位置。

2）在边坡的施工过程中，各岩石力学参数随之发生动态变化，其中弹性模量E变化最为显著，黏聚力C次之，参数上的变化能反映边坡变形的内在物理影响机制。

3）采用BP 神经网络耦合有限元模型对参数反分析，逼近实际监测结果，参数计算结果能有效反映实际情况。

4）通过计算位移值对比分析施工现场位移信息，揭示了岩体参数在施工过程中的动态变化过程，获得了边坡稳定与变形的关系，该方法适用于深挖碎裂岩质边坡的稳定性分析。

5）采用施工与数值模拟结合，计算分析获得最优模型计算参数来进行稳定性分析，提供了一种有效可行的信息化施工方法。