基于物质点强度折减法的二维...质c-φ边坡稳定性分析图表_张鹏.pdf

书书书Journal of Engineering Geology工程地质学报10049665/2023/31(2)-0596-11张鹏，张绍和，刘磊磊 2023 基于物质点强度折减法的二维均质 c 边坡稳定性分析图表J 工程地质学报，31(2):596606 doi:1013544/jcnkijeg2020660Zhang Peng，Zhang Shaohe，Liu Leilei 2023 Two-dimensional homogeneous c slope stability analysis charts based on material point strength reductionmethodJ Journal of Engineering Geology，31(2):596606 doi:1013544/jcnkijeg2020660基于物质点强度折减法的二维均质 c 边坡稳定性分析图表*张鹏张绍和刘磊磊(有色金属成矿预测与地质环境监测教育部重点实验室(中南大学)，长沙 410083，中国)(中南大学地球科学与信息物理学院，长沙 410083，中国)摘要安全系数是判断边坡稳定性最常用的指标。应用传统极限平衡法求解边坡安全系数时，没有考虑土体内部的应力应变关系，不能准确、真实地描述边坡土体破坏时的应力场、位移场及塑性区的发展情况。而应用有限元法虽然可以解决上述问题，但在计算过程中，会因网格畸变导致计算不收敛，造成求解困难。为此，本文基于物质点强度折减法，建立了安全系数无量纲参数 F/tan 与土体无量纲参数 之间的关系，并绘制了求解均质边坡安全系数的稳定性图表。该方法无需任何迭代过程即可快速、方便地确定边坡安全系数。以文献中常见的 4 个典型边坡算例为例，将本文计算结果与文献中其他边坡稳定性图表进行对比，验证本文方法的准确性。结果表明，用本文研制的边坡安全系数分析图表计算获得的边坡安全系数与其他方法获得的安全系数计算结果相接近，吻合较好，证明该方法是一种有效且保守的估计边坡安全系数的计算方法，为快速评估边坡安全系数提供了一种新的途径。关键词边坡稳定性分析;稳定性图表;物质点法;强度折减法中图分类号:TU42文献标识码:Adoi:1013544/jcnkijeg2020660*收稿日期:20201228;修回日期:20210506基金项目:国家自然科学基金项目(资助号:41902291)，湖南省自然科学基金项目(资助号:2020JJ5704)This research is supported by the National Natural Science Foundation of China(Grant No 41902291)and the Natural Science Foundation of HunanProvince(Grant No 2020JJ5704)第一作者简介:张鹏(1997)，男，硕士生，主要从事滑坡过程大变形模拟研究 E-mail:pengzhangcsueducn通讯作者简介:刘磊磊(1987)，男，博士，副教授，主要从事地质灾害防治与风险控制研究 E-mail:csulllfoxmailcomTWO-DIMENSIONAL HOMOGENEOUS c SLOPE STABILITY ANALYSISCHATS BASED ON MATEIAL POINT STENGTH EDUCTION METHODZHANG PengZHANG ShaoheLIU Leilei(Key Laboratory of Metallogenic Prediction of Nonferrous Metals and Geological Environment Monitoring，Ministry of Education，School of Geosciences and Info-Physics，Central South University，Changsha 410083，China)(School of Geosciences and Info-Physics，Central South University，Changsha 410083，China)AbstractFactor of safety(FOS)is an important parameter to evaluate slope stability Two common methods usedfor evaluating FOS are the limit equilibrium method(LEM)and the finite element method(FEM)Since the stress-strain relationship is not considered，the LEM cannot truly reflect the development of stress field，displacement fieldand plastic zone during soil failure On the other hand，although the constitute model is considered in FEM，it maysuffer from nonconvergence problem due to grid distortion in the calculation process To solve this problem，thestrength reduction material point method(SMPM)is proposed to evaluate the slope stability The relationship be-tween dimensionless soil parameters and FOS is subsequently constructed and plotted by different curves to obtainthe stability analysis chart This method can quickly and conveniently determine the FOS of slope without any itera-tive process The reliability of the method is verified by four typical examples The results show that the FOS calcu-lated by the proposed chart is close to that calculated by other methods The proposed method is conservative andeffective.Key wordsSlope stability analysis;Stability chart;Material point method;Strength reduction method0引言滑坡是自然界中最常见的地质灾害(林松等，2019)，这种地质灾害发生突然，难以提前准确预测，给滑坡滑动路径中及滑坡下部人们的生产生活带来很大的安全隐患，常造成严重的人员伤亡和财产损失(Zhang et al，2014)。在评估边坡稳定性时，安全系数是判断边坡稳定性的常用指标(刘文红等，2016;苏永华等，2019)。一般而言，使用有限元法和极限平衡法等方法求解边坡安全系数时，首先需要在边坡稳定性分析软件或程序代码中进行必要的前处理(如建立边坡模型，输入土体参数，根据边坡状态设置边界条件，施加重力或均布荷载及划分网格等);前处理完成后，才能进行稳定性计算。该方法虽然应用比较成熟，且能在短时间内获得边坡稳定安全系数，但针对同类型的边坡，依然需要重复上述前、后处理过程，因而分析过程有待进一步优化。稳定性图表作为快速分析和评价边坡稳定性的有效工具，则可避免上述必要的前处理过程，在很大程度上减少了边坡稳定性分析的计算成本，因此具有一定的实际工程意义。1937 年 Taylor(1937)首次提出稳定性图表这一概念。随后，Baker(2003)在 Taylor 公式基础上定义了临界滑移面，建立了稳定数与临界滑移面之间的关系，并以简单图表表示。Steward et al(2011)回顾了 Taylor 提出的两种边坡稳定性设计图表，并分析绘制了具有不同几何形状和土性参数的各种稳定性图表。Sahoo et al(2019)在 Majumdar(1971)的基础上绘制了考虑水平和竖向地震力影响的泰勒稳定性图表。Vo et al(2017)在 Morgenstern et al(1960)工作的基础上，利用滑移线理论对由非饱和土和非均质土组成的曲线斜坡进行稳定性分析，并给出了非均质非饱和土中曲线边坡的稳定性分析图表。随着时代的进步和研究的深入，边坡稳定性分析理论也在不断地发展和完善，领域内的许多专家学者在传统极限平衡法和极限分析法的基础上提出了一系列边坡稳定性分析图表。在极限分析框架下，He et al(2020)结合修正的屈服准则，对边坡安全系数进行解析求解，得到了边坡安全系数的上限解，并给出了边坡在不同地震加速度作用下的稳定性图表。Michalowski(2002)基于极限分析运动学方法，给出了孔隙水压力及地震力作用下黏性土坡稳定性分析图表。为了避免安全系数陷入极小值，Tang et al(2015)应用强度折减法和上限分析定理，建立了简单均质边坡安全系数的目标函数，得到了严格的安全系数上限，并将线性序列二次规划法(SQP)与随机游走法相结合，对均质各项同性边坡安全系数的目标函数进行优化，绘制了简单均质边坡安全系数稳定性分析图表。但是，用传统方法求解边坡安全系数时，没有考虑土体的应力应变关系，仅关心岩土体最终整体滑动时的状态，这种计算方法不能准确、真实地再现土体失稳破坏时的应力场、位移场及塑性区的发展情况。针对传统边坡稳定性分析方法的不足，孙超伟等(2018)利用 UFIELD 平台对有限元分析软件ABAQUS 进行开发升级，并将场变量和有限元强度折减法相结合，建立了一套求解均质边坡安全系数的稳定性分析图表。年廷凯等(2012)将库水位变化及渗流条件加入 ABAQUS 计算程序中，建立了一体化的安全系数稳定性分析图表。与其他边坡稳定性分析方法相比，该方法能够很好地再现边坡破坏过程并充分考虑了土体内部的应力应变关系，不需要根据钻孔数据和室内试验人为假定边坡临界滑动面的形状和位置，便可求得任意类型边坡的临界滑动面和边坡处于极限状态的最小安全系数，降低了因滑动面漏划、错划等人为因素造成的误差。然而，这种数值方法在计算过程中会因网格畸变导致雅可比矩阵为负，使计算不收敛，造成求解困难。近年来发展的物质点法(MPM)采用欧拉网格和拉格朗日质点双重描述，通过背景网格与物质点间的交替映射完成数值计算，每个计算时间步均采79531(2)张鹏等:基于物质点强度折减法的二维均质 c 边坡稳定性分析图表用未变形的背景网格集成动量方程并进行空间导数的求解(孙玉进等，2015)，避免了因网格畸变(Lianet al，2013)而带来的数值求解困难。这种数值分析方法弥补了有限元法的不足，可以很好地解决边坡大变形问题，因此，越来越多的人开始将物质点法延伸到边坡领域，用于求解边坡稳定性问题(史卜涛等，2016;王安礼等，2016;王双等，2016)。但物质点方法计算效率低下，且目前没有成熟的商业软件，导致其应用具有一定的局限性。因此，为了更好地推动物质点法在边坡稳定性分析方面的发展，很有必要进一步建立边坡稳定性分析的物质点方法设计图表，进而为边坡稳定性快速分析评价提供另一种可能。为此，本文以典型二维均质 c 边坡为例，提出将物质点法与强度折减法相结合，并参考常见的参数取值范围，旨在建立安全系数无量纲参数F/tan 与土体无量纲参数 之间的关系，提出求解二维均质 c 边坡安全系数的稳定性分析图表，以便用物质点法进行边坡稳定性分析时，可快速、方便地获得边坡安全系数。最后，通过 4 个典型实例验证了该方法的可靠性。1物质点强度折减法边坡稳定性分析基本原理与有限元强度折减法基本原理类似，物质点强度折减法亦是采用不断减小边坡土体抗剪强度参数的方式，迫使边坡达到失稳破坏状态，从而利用物质点法数值计算程序自动获取边坡滑动破坏面和安全系数(赵尚毅等，2002;Zhao et al，2012)。该方法具有有限元法的各种优势，包括考虑边坡土体的应力应变关系和不需要人为假定滑动面的形状和位置，但不存在网格畸变使计算不收敛等问题。其基本原理主要包括物质点法、强度折减理论和边坡失稳判据，分别介绍如下。1.1物质点法1994 年，Sulsky et al(1994)为了将流体力学中的质点网格法(FLIP)应用到固体力学领域，提出了物质点法。物质点法采用拉格朗日质点和欧拉网格双重描述，将材料区域离散成一系列质点(如图 1所示)(史卜涛等，2016)，每个质点代表一块区域并携带该区域的所有物质信息，而网格在空间固定，不携带任何信息，仅用于动量方程的求解和空间导数的计算(张雄等，2011)。图 1物质点法示意图(张雄等，2013)Fig 1Schematic diagram of MPM:(a)material region，(b)particle representative region and(c)discretization of thematerial points(Zhang et al，2013)a 材料区域;b 质点代表区域;c 物质点离散基于更新的拉格朗日控制格式，其动量方程为:ijxj+bi=ui(1)式中:为当前材料质量密度;xj为空间坐标;ij为Cauchy 应力;bi为单位质量的体积力;ui为位移。其中:边界条件为:(njiji)=tiviu=vi(2)式中:S和Su分别为给定面力边界条件和给定位移边界条件;nj为边界 S的外法线单位向量;ti为作用在边界S的面力;vi为位移边界Su的速度。更新拉格朗日格式等效积分弱形式(也称为虚功方程)为:uuidV+sijui，jdV biuidVtsiuidV=0(3)式中:V 为现实构型的体积;sij为比应力;?tsi为比边界面力;ui为虚位移。物质点法计算可分为 3 个阶段，即初始化阶段、更新拉格朗日计算阶段和映射阶段。在初始化阶段，如图 2a 所示，区域中的所有物质点将自身携带的信息传递到背景网格节点上，计算背景网格节点质量、动量等信息，为下阶段做准备。更新拉格朗日计算阶段(图 2b)则在背景网格节点上积分动量方程，得到新的速度、加速度、位移等信息。最后，在映射阶段(图 2c)，基于上阶段计算结果，更新所有物质点信息，该时间步计算完成，网格恢复原状，进行下一时间步计算(图 2d)。需要指出的是，边坡稳定性分析属于静力问题，而物质点法为动力分析，在利895Journal of Engineering Geology工程地质学报2023图 2物质点法计算流程Fig 2Calculation process of MPM用物质点法进行边坡稳定性分析时，需要考虑动力分析引起的数值振荡问题。为此，本文在进行物质点边坡稳定性分析计算前，通过分步加载重力加速度，采用弹塑性模型模拟边坡初始应力场，以消除过大应力引起的数值振荡(史卜涛等，2016)。1.2强度折减理论强度折减法的基本原理是保持土体的重力加速度为常数，将边坡土体强度参数除以相同的折减系数 K，得到当前折减系数下的黏聚力 c1和内摩擦角1，并把折减后的土体参数作为新的计算参数代入物质点法程序中进行计算。经过数次折减，边坡土体达到极限破坏状态失稳(韩伟歌等，2019;刘康琦等，2020)。此时的折减系数定义为边坡的安全系数(Zienkiewicz et al，1975;赵尚毅等，2002;Zhao et al，2012)，相应的破坏面即为边坡失稳滑面。强度参数折减公式为(霍沿东等，2019):c1=c0K(4)1=arctantan 0K(5)式中:K 为折减系数;c0和 c1分别为折减前、后的黏聚力;0和 1分别为折减前、后的内摩擦角。1.3边坡失稳判据失稳判据的选择是求解边坡安全系数过程中的必要环节。在岩土工程实践中，使用最频繁的边坡失稳判据主要有下列 4 类(唐宇峰等，2016;涂义亮等，2018):(1)有限元计算不收敛(Ugai et al，1995;Griffiths et al，1999;郑颖人等，2004);(2)广义剪应变贯通(连镇营等，2001;刘鑫等，2019);(3)特征点位移突变(Zienkiewicz et al，1975);(4)形成贯通的塑性区。但工程界对有限元强度折减进行边坡稳定性分析的失稳判据并没有统一标准，采用的比较多的是静力计算不收敛准则。然而，因物质点法采用中心差分法显示求解(张雄等，2013)，满足 CFL 条件，数值计算是稳定的。根据 LAX 等价原理，这种差分格式是收敛的，因此不能使用数值计算不收敛作为边坡失稳判据。文献(史卜涛等，2016)指出，对于同一边坡，物质点法比有限元法形成的塑性区大，且贯通后，在很长一段时间内，边坡不会发生失稳破坏，因此本文使用坡顶点竖直方向的位移是否突变作为失稳判据。2物质点强度折减法边坡稳定性分析实施步骤为了理清本文物质点强度折减法进行边坡稳定性分析的思路和促进理解，本节详细介绍物质点强度折减边坡稳定性分析流程，具体如图 3 所示。从图 3 可以看出，物质点强度折减法边坡稳定性分析主要包括前处理(即边坡离散)、物质点强度折减法计算(即初始应力场生成、强度折减、物质点法计算)和后处理(即安全系数求解)等 3 个过程，并具体由以下步骤组成:步骤 1:前处理，即边坡域离散。因物质点法是将连续体离散成一系列质点，每个质点代表一块区域并携带了该区域包括质量、密度、应力、应变在内的所有物质信息，因此，应用物质点法进行边坡稳定性计算时，首先需对边坡进行离散。本文使用有限元分 析 软 件 ABAQUS 对 边 坡 进 行 建 模，利 用ABAQUS 集成的 python 2.7.0 和 MATLAB 2019b 获取边坡离散信息，具体包括:(1)在 ABAQUS 中建立边坡模型，划分网格，导出含有单元编号和网格节点坐标的边坡模型“inp”文件;(2)利用 ABAQUS 自带的 python 接口、数据处理软件 MATLAB 2019b 和导出的“inp”文件，将建立的边坡模型离散成物质点，并导出物质点坐标、面积和质量。步骤 2:生成边坡初始应力场。如前所述，物质点是动力分析方法，边坡稳定性分析是静力问题，用动力物质点法分析边坡静力问题时会出现数值振荡、计算不准确等问题(孙玉进等，2015)，因此，在进行边坡稳定性分析之前需要确定边坡的初始应力场(史卜涛等，2016)。本文采用弹性模型确定边坡初始应力场，且重力需分步加载(Andersen et al，99531(2)张鹏等:基于物质点强度折减法的二维均质 c 边坡稳定性分析图表图 3物质点强度折减法边坡稳定性性分析流程图Fig 3Flowchart for material point methodstrength reduction method2010)以消除过大应力引起的数值振荡。对于重力加速度(张巍等，2017)的加载如图 4 所示，02is按照 线 性 增 加 的 方 式 加 载 重 力 加 速 度，直 至9.8im s2;2is 后重力加速度不变。0 3is 整个加载过程土体采用线弹性模型，3is 后，土体采用服从Drucker-Prager(DP)准则的弹塑性本构模型。图 4重力加载方式(史卜涛等，2016)Fig 4Incremental gravitation(Shi et al，2016)步骤 3:物质点法计算。首先利用第 1.2 节所述的强度折减方法，按当前的折减系数对土体强度参数依次折减，获得当前折减系数下的土体强度参数 c1和 1，其次在步骤 2 生成的边坡初始应力场的基础上，将 c1和 1作为新的边坡土体强度参数，重新导入物质点法数值计算程序中进行边坡稳定性计算，并记录边坡特征点最大位移。步骤4:判断是否满足计算终止条件，否则继续折减边坡土体强度参数，重复步骤3，直至边坡失稳破坏。步骤 5:利用 1.3 节边坡失稳判据方法，绘制位移折减系数曲线图，获取边坡安全系数。3本文方法可行性分析为了说明本文方法的有效性，采用文献(史卜涛等，2016)中的算例对本文方法进行验证。该算例边坡的几何外形如图 5 所示，其中坡高 H=10 m，坡脚=45，膨胀角=0，土体重度=20ikNm3;土体黏聚力和内摩擦角分别为 c=12.38ikPa和=20，弹性模量 E=100iMPa，泊松比=0.30。采用 DP 屈服准则的弹塑性本构模型。在本次计算中，该边坡模型底部在 3 个方向均固定，用以模拟下部基岩;两侧、前面和后面都采用一端固定，一端自由的边界条件，用以模拟边坡土体在平面应变情况下的各种状态。图 5算例边坡几何尺寸Fig 5Slope geometry of example用本文物质点强度折减法计算该边坡安全系数，并与其他方法的结果进行对比，结果列于表 1。由表 1 可知，用极限平衡法计算出的该边坡安全系数为 1.00，用有限元强度折减法以塑性区贯通作为失稳标准计算出的安全系数为 0.98，以特征点位移突变作为失稳判据计算出的安全系数为 0.99，以计算不收敛为失稳判据得出的安全系数为 1.06，因此该边坡安全系数在 0.98 1.06 之间。而用本文物质点强度折减法计算出的该边坡安全系数为 1.00，在极限平衡法和有限元强度折减法计算的安全系数范围内，且与其他方法的计算结果吻合较好，从而验证了本文方法在计算边坡安全系数方面的准确性。为了进一步验证用本文物质点强度折减法计算边坡安全系数的准确性，绘制了不同折减系数下边006Journal of Engineering Geology工程地质学报2023表 1不同方法获得的边坡安全系数Table 1FOS obtained by different methods计算方法安全系数Spencer 极限平衡法(史卜涛等，2016)100SFEM(塑性区贯通)(史卜涛等，2016)098SFEM(特征点位移突变)(史卜涛等，2016)099SFEM(计算是否收敛)(史卜涛等，2016)106物质点强度折减法(特征点位移突变)(史卜涛等，2016)100本文方法100图 6不同折减系数下边坡特征点最大位移随迭代时步变化曲线Fig 6The curve of maximum displacement of slopecharacteristic points with iterative time step underdifferent reduction coefficients坡特征部位(坡顶点)最大位移随迭代时步变化曲线，如图 6 所示。由图 6 可知，当折减系数小于 1.00 时，特征点最大位移随迭代时间步的增长变化缓慢并趋向稳定。当折减系数大于 1.00 时，特征点最大位移随迭代时间步持续增长且增长速度越来越快，无趋向稳定的趋势。由此可知，当折减系数等于或小于 1.00时，边坡是稳定的，仅发生很小的位移，当折减系数大于 1.00 时，边坡已经失稳，发生了大位移。因此，折减系数 1.00 是边坡安全系数的界限值，可认为该边坡的安全系数为 1.00，该计算结果与其他方法获得的结果基本一致，进一步验证了本文方法计算边坡安全系数的准确性，证明该方法可用于边坡安全系数求解。为了进一步说明本文方法确定的边坡失稳滑面的有效性，图 7a 和图 7b 分别给出了同一折减系数下，用两种强度折减法计算出的该边坡等效塑性应变云图。在本次计算中，选用平面应变 8 节点的四图 7边坡塑性区分布Fig 7Plastic region of the slopea 有限元法塑性区分布;b 物质点法塑性区分布边形网格单元在 ABAQUS 软件中进行有限元计算，其中:土体采用 Mohr-Coulomb(MC)强度准则;物质点法计算采用 8 节点立方体单元作为背景网格，单元尺寸为 1 m1 m1 m。需要说明的是，物质点间距会影响边坡稳定性计算的时间和效率，因此，折中考虑，多次试验后本文最终选用的物质点间距为0.5im，即每个背景网格单元采用 8 个物质点，共得到 770 个物质点。图 7b 是用物质点法获得的边坡有效塑性应变分布图，由该图可知，该边坡在当前折减系数下已经失稳，且有向右侧滑动的趋势，而左侧设置了一端固定，一端自由的对称边界条件，限制了其水平滑动，所以边坡顶部产生了变形;坡脚处由于边界条件突变，产生了较大的塑性变形且塑性区贯通至坡顶。由图 7a 和图 7b 可知，用两种强度折减法得到的该边坡塑性区范围基本一致，说明用本文方法计算得到的塑性应变及塑性区范围合理有效。综上所述，本文提出的物质点强度折减法的塑性区分布范围、塑性应变及安全系数与其他方法计算出的结果基本一致，说明该方法的计算结果准确可靠，可用于边坡安全系数求解，也为后续均质边坡稳定性图表设计提供了保障。10631(2)张鹏等:基于物质点强度折减法的二维均质 c 边坡稳定性分析图表图 8边坡二维计算模型Fig 8Numerical 2D slope model4二维均质 c 边坡稳定性分析图表本节采用所提出的物质点强度折减法研制二维均质 c 边坡稳定性分析图表。二维均质边坡稳定性分析涉及到的几何参数有坡高 H、坡脚，材料参数有重度、黏聚力 c、内摩擦角、弹性模量 E 和泊松比(高冯等，2020)。如果单一考虑各个因素对边坡稳定性的影响，在数值计算及后期图表处理存在较大困难。因此，为了方便图表表达，将上述土体参数分类整合并做无量纲化处理，得到了安全系数无量纲参数F/tan 与土体材料无量纲参数=c/Htan。这种表达式的优点在于参数 与安全系数 F 相对独立，且涵盖坡体的基本土性参数，大大简化了图表参数表达式(年廷凯等，2012)。考虑的二维边坡数值计算模型如 图 8 所示。由于边界效应的影响，模型坡顶到左边界的距离L1、坡脚到右边界的距离 L2均为坡高 H 的 2.0 倍，下边界距离坡脚的垂直高度 D 为坡高 H 的 1.0 倍。参考文献中关于边界约束条件的选取原则，在本次稳定性分析图表计算中，该计算模型底部在 x、y 和 z方向上均使用固定边界条件，以模拟边坡下部基岩;边坡两侧、前面和后面均设置一端固定、一端自由的对称边界条件，用以模拟平面应变情况下边坡土体的各种状态。计算中土体采用服从 DP 屈服准则的弹塑性本构模型，采用六面体 8 节点单元的背景网格，单元尺寸为 1 m1im1im，折中考虑计算效率与精度物的影响，物质点间距与前述算例一样，取为0.5im;计算时间步长为 0.465103s，总计算时长为 20is。参考常见的参数取值范围，本文设计的坡高 H 分别为 5im、10im、15im、20im;坡角分别为15、30、45、60、75、90;土体重度 分别为16ikN m3、18ikN m3、20ikN m3;黏聚力 c 分别为5ikPa、10ikPa、15ikPa、20ikPa、25ikPa、30ikPa、35ikPa、40ikPa;内摩擦角 分别为 5、10、15、20、25、30、35、40。张培文等(2006)指出泊松比 的取值对边坡计算形成的塑性区范围有一定影响，且计算结果显示，的取值与塑性区范围呈反比关系，即 越大，塑性区范围越小;但数值计算结果显示，=0.1 和=0.49 计算得到的安全系数是相同的，说明在安全系数计算过程中，可以忽略 对计算结果的影响;赵尚毅等(2005)也指出，弹性模量E 对边坡安全系数的计算结果和形成的塑性区范围只产生很小的影响，且对安全系数计算结果的影响不显著。在本次边坡稳定性分析图表计算中，采用特征点位移突变判据，而没有使用塑性区贯通边坡失稳判据，因此，在本次二维均质 c 边坡稳定性分析图表计算过程中，视变形指标弹性模量 E 和泊松比 为常数，即弹性模量 E=100iMPa，泊松比=03。基于上述参数，在初始应力场的基础上，用物质点强度折减法对按照上述边坡模型及计算参数设计的边坡分别进行稳定性计算，并将计算结果线性拟合，绘制成F/tan-关系曲线，如图 9 中实线所示。为了验证本文研制的图表的有效性，图 9 也给出了文献中的相关稳定性分析图表，并用虚线表示。其中:图 9a 和图 9b 中文献结果采用有限元强度折减法进行边坡安全系数计算，图 9c图 9e 中文献则采用极限分析法获得边坡安全系数。由图 9 可以看出:本文获得的稳定性图表与文献(Michalowski，2002;Michalowski，2010;Gao etal，2013)中的稳定性图表在数值和趋势上基本一致，而与文献(孙超伟等，2018;Sun et al，2019)中的稳定性图表在趋势上基本一致，但在数值上差别较大。具体而言，由图 9a 和图 9b 可知，当边坡坡脚 为 15和 30时，本文获得的安全系数曲线与文献(孙超伟等，2018;Sun et al，2019)中获得安全系数曲线在数值和趋势上基本一致，而当边坡坡脚 为 45、60、75和 90时，本文获得的安全系数曲线与文献(孙超伟等，2018;Sun et al，2019)中获得安全系数曲线虽在趋势上基本一致，但在数值上差异较大。主要原因是图 9a 中文献(Sun et al，2019)采用有限元数值计算程序迭代不收敛作为边坡失稳判据求解边坡安全系数，图 9b 中文献(孙超伟等，2018)采用有限元数值程序计算过程中形成的塑性区从坡脚贯通至坡顶作为边坡失稳判据求解边坡安全系数，本文物质点强度折减法则是以特征点位移突变作为边坡失稳判据来求解边坡安全系数。而采用不同的失稳判据，所得到的安全系数大小会存在明显差异(史卜涛等，2016)。同时，计算方法本身的不同，亦会导致边坡稳定性分析结果的差别。因206Journal of Engineering Geology工程地质学报2023图 9与其他文献对比图Fig 9Comparison with other literatures此，本文获得的安全系数曲线与图 9a 和图 9b 对比文献中的安全系数曲线在数值上存在一定差异。然而，因边坡安全系数受几何参数和材料参数共同控制，且稳定性图表是通过安全系数无量纲参数与边坡内摩擦角的乘积来获得边坡安全系数，这种差异在安全系数结果上并不大。总之，上述不同图表之间的一致性间接说明了本文研制的基于物质点强度折减法的边坡稳定性分析图表的有效性。5算例验证为直接验证所提图表的有效性，本节从文献中选取了 4 个不同算例，对比了物质点强度折减法得到的安全系数与本文图表获取的安全系数。同时，作为进一步验证，文献中的方法也一并给出进行参考。5.1算例 1算例 1 为一均质土坡，坡高 10im，坡脚=45，土体重度=20ikN m3，黏聚力 c=12.38ikPa，内摩擦角=20，c/Htan=0.17。用本文获得的安全系数图表计算该边坡安全系数并与其他方法的计算结果作对比，计算结果如表 2 所示。5.2算例 2算例 2 为一均质土坡，坡高 8im，坡脚=30，土体重度=20ikNm3，黏聚力 c=20ikPa，内摩擦30631(2)张鹏等:基于物质点强度折减法的二维均质 c 边坡稳定性分析图表表 2各种方法得到的安全系数Table 2Values of factor of safety obtained by various methods计算方法安全系数Spencer 极限平衡法(史卜涛等，2016)100有限元强度折减法(史卜涛等，2016)098图表法(Tang et al，2015)095图表法(Michalowski，2002)100图表法(Michalowsk et al，2011)098物质点强度折减法098图表法(本文)097角=20，c/Htan=0.22。利用本文建立的安全系数图表计算该边坡安全系数，并与其他方法的计算结果作对比，计算结果如表 3 所示。表 3各种方法得到的安全系数Table 3Values of factor of safety obtained by various methods计算方法安全系数图表法(Tang et al，2017)185图表法(Tang et al，2015)206图表法(Michalowski，2002)233图表法(Michalowski，2010)237图表法(孙超伟等，2018)231极限平衡法186物质点强度折减法195图表法(本文)1985.3算例 3算例 3 为一均质土坡，坡高 8im，坡脚=45，土体重度=18.5ikNm3，黏聚力 c=20ikPa，内摩擦角=15，c/Htan=0.50。用本文获得的安全系数图表计算该边坡安全系数并与其他方法的计算结果作对比，计算结果如表 4 所示。表 4各种方法得到的安全系数Table 4Values of factor of safety obtained by various methods计算方法安全系数图表法(Michalowski，2010)134图表法(Tang et al，2015)134图表法(Michalowski，2002)131图表法(Tang et al，2017)125图表法(孙超伟等，2018)134物质点强度折减法128极限平衡法135图表法(本文)1255.4算例 4算例 4 为一均质土坡，坡高 10im，坡脚=30，土体重度=17ikNm3，黏聚力 c=10ikPa，内摩擦角=20，c/Htan=0.162。用本文获得的安全系数图表计算该边坡安全系数并与其他方法的计算结果作对比，如表 5 所示。表 5各种方法得到的安全系数Table 5Values of factor of safety obtained by various methods计算方法安全系数图表法(Tang et al，2015)140图表法(Tang et al，2017)118图表法(Michalowski，2002)131图表法(Michalowski，2010)131图表法(孙超伟等，2018)129极限平衡法131物质点强度折减法119图表法(本文)120由表 2至表 5 可知，本文研制的二维均质边坡稳定性分析图表计算出的边坡安全系数与其他图表法和极限平衡法计算出的结果基本一致，略低于其他方法的计算结果，但总体较为接近。因本文提出的稳定性分析图表是在大量数据的基础上线性拟合而成，且安全系数是通过安全系数无量纲参数与边坡内摩擦角正切值的乘积获得的，所以，直接用物质点强度折减法计算出的边坡安全系数与用本文提出的安全系数图表法计算出的安全系数不一致，但这种结果上的差别在误差允许范围内。综上所述:本文研建的二维均质 c 边坡稳定性分析图表计算出的安全系数与其他图表法和极限平衡法计算出的结果在数值上基本一致，吻合较好，证明了本文边坡稳定性分析图表的准确性。6结论本文在物质点强度折减法的基础上，研建了一套求解二维均质 c 边坡安全系数的稳定性分析图表，并详细介绍了物质点强度折减法基本原理和实施流程，得到以下结论:(1)基于物质点法，在初始应力场的基础上，将物质点法与强度折减理论相结合，提供了一种求解边坡安全系数的有效途径。具体计算过程为:边坡域离散，分步加载重力以消除突加荷载产生的数值406Journal of Engineering Geology工程地质学报2023振荡，确定边坡初始应力场;折减强度参数;用物质点法计算边坡特征点位移;根据计算结果求解边坡安全系数。(2)用本文提出的边坡稳定性计算方法对一个典型边坡算例进行稳定性计算，结