边坡稳定性分析作为岩土工程的重要内容，其核心内容主要包括计算边坡安全系数和确定边坡潜在滑动面^[1]。目前分析边坡稳定性的数值方法主要包括有限元法、离散元法、边界元法、物质点法、数值流行法、光滑粒子流体动力学法和连续-离散耦合法等^[2-9]。

数值方法中有限元法更为成熟，其在边坡工程中应用较为广泛，不少学者^[10-11]采用有限元法模拟实际边坡工程失稳破坏过程，建立了边坡安全系数计算方法。赵凌峰等^[12]利用有限元强度折减法计算在湿干-冻融循环条件下渠道的稳定性，根据其滑动破坏的原因进行针对性改造。彭普等^[13]基于强度折减法对某经典边坡案例进行了系统分析，提出一种边坡可靠度分析的新方法。由于边坡岩土体是一种复杂的非连续介质体，采用有限元难以准确描述其细观结构的变化，而离散元法作为描述非连续介质的方法，可以很好地从细观角度上模拟边坡的破坏过程^[14]。周健等^[15]将强度折减法的理念引入颗粒流离散单元法中，求解安全系数时相较于传统的极限平衡法不需要假定滑移面的位置和形状。金磊等^[16]引入离散元强度折减法并提出基于能量演化的方法来判断三维土石混合体边坡的失稳破坏条件。此外，边界元法、物质点法、数值流行法、光滑粒子流体动力学法等在研究边坡破坏过程和稳定性方面也得到了一些应用，对描述边坡破坏机制和评价边坡稳定性具有指导作用^[6-9]。

离散元法受限于颗粒数量与计算效率的缺点，往往不得不对颗粒粒径进行放大^[15]或对边坡模型尺寸进行缩小^[16]，因而难以完全还原实际边坡工程情况。连续-离散耦合方法选取重点关注区域为离散域，其余区域为连续域，通过数据交换通道进行模拟计算，极大地减少了颗粒数量，可以很好地弥补离散元计算效率慢的缺点。张铎等^[17]基于连续-离散耦合方法分析某尾矿坝边坡冲填前后的宏、细观力学特征，验证了耦合方法的可行性。CHEN等^[18]通过对比连续-离散耦合方法和其他已知方法的结果，证明了耦合方法在模拟重力和地震激励作用下非粘性边坡破坏的适用性。ANTOLINI等^[19]指出连续-离散耦合方法能够真实地模拟边坡破坏的不同阶段，从而用于边坡风险评估和管理。连续-离散耦合方法发展至今，大部分学者利用此方法进行边坡破坏模式和机理等定性研究，而对安全系数和滑动面方面的定量研究颇为少见。严琼等^[20]提出强度折减理论可用于基于连续-离散耦合方法的数值模型中，但并未给出具体的实施步骤以及安全系数的计算。因此，为了定量分析边坡的稳定性，开发一种适用于连续-离散耦合边坡模型的强度折减法至关重要。

强度折减法其关键问题之一是选用适合的边坡失稳判据，传统的失稳判据主要分为3种^[21]：数值计算不收敛^[21]、特征点位移突变^[22]、塑性区贯通^[23]。为了解决环境复杂多变的边坡工程问题，多名学者^[24-25]提出了若干新的边坡失稳判据，比如能量突变失稳判据。此外，涂义亮等^[26]证明了传统3种失稳判据在理论上具有统一性。目前已有的失稳判据多用于有限元强度折减法中，其在连续-离散耦合方法中的适用性尚且不知，探寻适用于连续-离散耦合强度折减法的失稳判据亦是开发该方法的一项重要内容。

已有研究表明^{[4, 26-27]}，边坡稳定性分析的精度与数值方法使用的合理性直接相关，其影响因素和提高方法较多。杨智勇等^[4]提出考虑边坡多失效模式的区域概率风险分析方法，更准确地量化边坡破坏的风险。涂义亮等^[26]研究表明有限元数值计算提高精度的根本办法在于兼顾计算效率的同时，严格收敛标准以及精细化网格。林珊等^[27]提出用于边坡稳定性分析的虚单元强度折减法，此方法的计算精度与单元形状无关，从而提高了安全系数的准确性。在连续-离散耦合边坡的研究中，离散域选取对边坡稳定性计算精度的影响尚不明确，目前这方面还缺少系统性地研究。

综上所述，为了解决基于连续-离散耦合方法的边坡稳定性定量分析问题，本文首先开发并验证一种连续-离散耦合强度折减法，其次从细观角度上分析边坡失稳破坏机理，然后探寻适用于耦合强度折减法的失稳判据，最后讨论离散域位置、大小及形状对边坡安全系数计算精度的影响，提出离散域的选取原则和方法。研究成果可为耦合强度折减法的精确应用提供一定参考。

1   连续-离散耦合强度折减法的开发

1.1   连续-离散耦合方法的基本原理

连续-离散耦合是指把一个已设置好的离散域嵌入到另一个连续域预留的空域中，并且在连续-离散耦合交界位置建立耦合分析数据交换通道，其基本原理如图1所示。在连续-离散耦合迭代计算过程中，连续域节点速度通过耦合边界节点传递至离散域，离散域接触力通过耦合边界节点传递至连续域，两者通过耦合边界不断传输交换力与速度的信息，其信息交换流程如图2所示。因此，该方法保证了离散域与连续域计算数据的一致性与连续性，进而实现分别从宏、细观角度上综合分析介质的力学运动过程^[20]。

图  1  连续-离散耦合方法的基本原理示意图

Figure  1.  Basic principle diagram of continuous-discrete coupling method

下载: 全尺寸图片 幻灯片

图  2  连续-离散耦合方法的信息交换流程

Figure  2.  Information exchange process of continuous-discrete coupling method

下载: 全尺寸图片 幻灯片

目前连续-离散耦合问题的核心主要有两点^[28]：① 标定离散域颗粒细观参数后其宏观表现要与连续模型的宏观性质一致；② 连续-离散耦合模型在耦合边界上要保证力和速度的连续性。

1.2   连续-离散耦合强度折减法的提出

为了解决基于连续-离散耦合方法的边坡稳定性定量分析问题，结合连续-离散耦合原理和强度折减法，现开发一种连续-离散耦合强度折减法，基本原理如图3所示，主要包括以下三个模块：

1)建立连续-离散耦合分析模型。耦合分析模型的建立方法总体上与常规数值模拟方法相似，主要在离散域设置、材料本构和参数等方面有所区别。边坡稳定性分析中，潜在滑动面、软弱夹层等位置应为重点关注区域，然后根据研究内容选取对应的重点关注区域为离散域，其他区域为连续域，在耦合交界处建立数据交换通道，进而从细观角度上分析边坡的稳定性。离散元颗粒细观参数直接决定了连续-离散耦合效果的优劣，应基于材料力学试验标定合适的离散元颗粒本构和细观参数。

2)建立宏、细观强度折减系数的关系。相比于有限差分和离散元强度折减法，耦合强度折减法更为复杂，必须同时对连续域和离散域的材料参数进行折减。但是，由于材料宏观参数与细观参数物理意义不一致，不存在一一对应关系，如何等效地折减宏、细观材料参数成为耦合强度折减法的关键。为此，提出了一种宏、细观强度折减系数关系的建立方法。该方法基于离散元三轴试验或其他力学试验测量细观强度参数折减${F_{{\rm{DEM}} - i}}$后的宏观强度参数${\tau _i}$，然后由式(1)反算出相匹配的宏观强度折减系数${F_{{\rm{FDM}} - i}}$，多次调整细观强度折减系数${F_{{\rm{DEM}} - i}}$，重复上述步骤，进而建立两者的关系式${F_{{\rm{FDM}}}} = f({F_{{\rm{DEM}}}})$，该关系式则可用于后续控制宏、细观强度参数的等效折减：

式中：${\tau _0}$为有限差分三轴试验在某围压下得到的抗剪强度；${\tau _i}$为离散元三轴试验在某围压下得到的抗剪强度。

3)进行连续-离散耦合强度折减计算。为计算边坡安全系数和滑动面，按照宏、细观强度折减系数的关系式，同时折减耦合模型的宏、细观强度参数，计算求解，逐步增大强度折减系数，重复上述步骤直至边坡达到失稳判据状态，从而输出边坡安全系数$F = {F_{{\rm{FDM}}}}$，获得宏、细观模拟结果。

图  3  连续-离散耦合强度折减法的基本原理图

Figure  3.  Basic theory diagram of the continuous-discrete coupling strength reduction method

下载: 全尺寸图片 幻灯片

2   连续-离散耦合强度折减法的验证

2.1   建立边坡的连续-离散耦合分析模型

基于提出的连续-离散耦合强度折减法，选取文献[26-29]中的经典边坡算例，采用FLAC^3D 6.0数值软件建立边坡连续-离散耦合模型，如图4所示。由图4可知，离散域设置在边坡滑动面附近，其大小为8 m×8 m，其余区域为连续域，边坡模型宽度为2 m，坡顶处设有位移监测点。建立边坡连续-离散耦合模型的主要步骤：首先，生成离散域颗粒试样，采用“落雨法”进行下落成样，对试样上方的墙体施加速度压缩试样至规定大小；其次，将试样四周墙体伺服至所在区域的初始地应力大小，伺服后离散域形状和大小发生一定的变化；最后，生成有限差分边坡模型进行初始地应力平衡，删除预留空域并导入离散域模型，设置耦合墙作为数据交换通道。

图  4  边坡模型尺寸　/m

Figure  4.  Size of slope model

下载: 全尺寸图片 幻灯片

建立有限差分和离散元三轴压缩试验模型，进行离散元细观参数的标定，如图5所示。其中，有限差分(FLAC^3D)网格采用边坡算例^[26-29]给定的Mohr-Coulomb本构模型，材料参数如表1所示；离散元(PFC^3D)颗粒间接触采用平行黏结模型。设置三轴压缩试验的围压为200 kPa，采用“试错法”不断改变离散元颗粒细观参数，最后得到相匹配的应力-应变曲线，如图5所示。结果表明离散元模型的应力-应变曲线峰值主应力与有限差分模型的接近，其应力变化趋势一致，具有较好的匹配度，标定的离散元颗粒细观参数如表2所示。同时，离散元模型的峰值强度所对应的轴向应变大于有限差分模型，造成这种现象的主要原因是受限于颗粒数量和计算效率，离散元颗粒半径不能过小，从而导致其轴向应变会因孔隙在加载过程中的调整而变大。这种误差将导致耦合边坡模型的变形略微大于有限差分边坡模型，且对边坡安全系数计算结果影响很小。

图  5  大主应力-应变曲线标定对比

Figure  5.  Comparison of large principal stress-strain curves

下载: 全尺寸图片 幻灯片

表  1  连续单元的宏观力学参数^{[26, 29]}

Table  1.  Macro-mechanical parameters of continuous elements

弹性模量 E/MPa	泊松比 ν	黏聚力c/kPa	内摩擦角 φ/(°)	重度 γ/(kN/m3)
100	0.3	42	17	20

下载: 导出CSV 

| 显示表格

表  2  离散颗粒的细观力学参数

Table  2.  Meso-mechanical parameters of discrete particles

颗粒半径 R/m	颗粒密度 ρ/(kg/m3)	有效模量 E*/MPa	刚度比 k*	颗粒间摩擦 系数μ	平行黏结有效 模量$\overline E {^ \text{*} }$/MPa	平行黏结 刚度比$\overline k ^*$	平行黏结法向 强度${\overline \sigma _{\rm{c}}}$/kPa	平行黏结 粘聚力$\overline c$/kPa	平行黏结 摩擦角$\overline \phi$/(°)
0.10~0.16	2000	8	1.8	0.32	32	1.8	95	90	25

下载: 导出CSV 

| 显示表格

2.2   建立边坡的宏、细观强度折减系数关系

根据提出的耦合强度折减法，设置三轴压缩试验围压为200 kPa，此时${\tau _0}$=171. kPa，细观强度折减系数${F_{{\rm{DEM}}}}$从1~40逐步增大，由式(1)计算得到对应的宏观强度折减系数${F_{{\rm{FDM}}}}$，如图6所示。通过拟合得到两者满足数学关系式(2)，拟合相关系数R²=0.99968。由式(2)可知，宏观强度折减系数${F_{{\rm{FDM}}}}$与细观强度折减系数${F_{{\rm{DEM}}}}$并非呈线性关系，而是呈指数关系。出现这种现象的主要原因是：随着摩擦系数折减程度的增大，其对颗粒间的摩擦力影响越来越小，最终导致其宏观强度逐渐降低并趋于稳定：

式中：A、B、C均为拟合参数，其中A=1.686 31，B=−0.682 93，C=0.3164；当${F_{{\rm{DEM}}}}$=1.0时，${F_{{\rm{FDM}}}}$=1.0。

最后，设置${F_{{\rm{FDM}}}}$为1.00~1.35，利用式(2)计算出${F_{{\rm{DEM}}}}$，同时折减耦合模型的宏、细观强度参数，计算直至边坡达到失稳判据状态，绘制坡顶竖向位移与宏观强度折减系数${F_{{\rm{FDM}}}}$关系的拟合曲线，并对曲线左右两端分别作两条切线，如图7所示。由图7可知，当${F_{{\rm{FDM}}}}$=1.2678时，坡顶位移开始急剧增加，故以特征点位移突变为判据^[22]，则边坡安全系数为1.2678。

图  6  宏、细观强度折减系数关系的拟合曲线

Figure  6.  Fitting curve of the macro strength reduction factor vers the meso-one

下载: 全尺寸图片 幻灯片

图  7  坡顶竖向位移与宏观强度折减系数关系的拟合曲线

Figure  7.  Fitting curves of the vertical displacement at the slope top vers the macro strength reduction factor

下载: 全尺寸图片 幻灯片

2.3   耦合计算位移连续性验证

连续-离散耦合的关键问题之一就是要在耦合边界处保证力和速度的连续性。图8为边坡耦合计算收敛的位移云图(安全系数1.2678)。由图8可知，离散域和连续域的最大、最小位移值接近，耦合边界上位移具有很好的连续性；耦合边界左上角附近出现明显的分层现象，虚线区域及其附近区域为边坡的潜在滑动面。因此，耦合计算位移的连续性证明了耦合强度折减法的可行性。

2.4   边坡安全系数的验证

为验证基于耦合强度折减法计算所得边坡安全系数的准确性，分别采用有限差分强度折减法、离散元强度折减法、Janbu法、简化Bishop法计算边坡安全系数，计算结果如表3所示。耦合强度折减法与其他安全系数计算方法计算结果极为接近，差距最大为5.39%，表明此方法计算的安全系数可靠，具有较高的计算精度。

图  8  边坡耦合计算收敛的位移云图(安全系数1.2678)

Figure  8.  Displacement contour in slope using coupling method (when safety factor is 1.2678)

下载: 全尺寸图片 幻灯片

表  3  各算法边坡安全系数对比

Table  3.  Comparison of slope safety factors of each algorithm

计算方法	耦合强度 折减法	有限差分 强度折减法	离散元强度 折减法	Janbu法	简化 Bishop法
边坡安全系数	1.2678	1.227	1.271	1.224	1.203
差距/(%)	—	3.33	0.25	3.58	5.39

下载: 导出CSV 

| 显示表格

2.5   潜在滑动面的验证

为确定连续-离散耦合方法在边坡工程中的适用性，有必要对比耦合强度折减法和有限差分强度折减法、离散元强度折减法、Janbu法、简化Bishop法计算的边坡潜在滑动面的分布情况，如图9所示。连续域为最大剪应变增量云图，离散域内虚线区域表示颗粒间的黏结已经破坏，剩余区域表示颗粒间黏结完好。由图9可知，在五种计算方法对应的安全系数下潜在滑动面均已贯通；耦合强度折减法的耦合边界处连续性很好，离散域内颗粒间黏结沿着潜在滑动面发生不同程度的破坏；相比于耦合强度折减法，Janbu法和简化Bishop法计算的潜在滑动面位置和形状与其高度一致；而有限差分强度折减法计算的潜在滑动面位置离坡面稍远，离散元强度折减法计算的潜在滑动面位置离坡面稍近，总体上耦合强度折减法与其它几种方法的结果基本一致。因此，上述结果证明了连续-离散耦合方法在边坡工程中具有较好的适用性。

图  9  不同安全系数计算方法下的边坡潜在滑动面示意图

Figure  9.  Diagram of potential slide face in slope using different calculation methods of safety factors

下载: 全尺寸图片 幻灯片

3   基于连续-离散耦合方法的边坡破坏机理和细观失稳判据分析

为研究边坡失稳破坏机理和细观失稳判据，结合前文所述的耦合方法，将边坡耦合模型计算至收敛，利用耦合强度折减法获取不同强度折减系数下离散域中颗粒接触力链和颗粒黏结状态的演化过程。

3.1   颗粒接触力链的演化和突变失稳判据

图10为不同强度折减系数下颗粒接触拉力链的演化过程，颗粒之间承受的径向拉力显示为拉力链，力链粗细表示拉力的大小。由图10可知，在初始状态A点时，离散域内的拉力链受力较为均匀，且方向趋于水平；随着${F_{{\rm{FDM}}}}$的增加，离散域内左上方的拉力链受力逐渐增大，出现了明显的分层现象；应力传递直至开始出现大面积破坏，右下方的拉力链与滑动面方向近乎垂直，边坡失稳破坏，此时右下方的拉力链也开始逐渐破坏。其中当${F_{{\rm{FDM}}}}$=1.235时，边坡处于临界破坏状态。上述结果说明边坡初始拉应力方向趋于水平，随着强度参数的不断劣化，边坡内部开始应力重分布，应力方向发生改变，各向异性程度增加，最终导致潜在滑动面贯通，边坡开始失稳破坏，滑动面向边坡内部延伸发展。C点为拉力链数量的突变点，其对应的${F_{{\rm{FDM}}}}$=1.235可作为边坡的安全系数，与前文以特征点位移突变为判据耦合计算的安全系数相比，差距仅为2.59%，故提出一种新的边坡失稳判据：拉力链数量突变判据。

图  10  不同强度折减系数下颗粒接触拉力链的演化过程

Figure  10.  Evolution of particle contact tension chain with different strength reduction factors

下载: 全尺寸图片 幻灯片

3.2   颗粒黏结状态的演化和突变失稳判据

图11为不同强度折减系数下颗粒黏结状态云图，图12为不同强度折减系数下颗粒黏结状态的演化过程。由图11和图12可知，随着${F_{{\rm{FDM}}}}$的增加，离散域内左上方的力链受力增大直至发生张拉、剪切破坏；黏结出现大面积破坏，左上方处的部分颗粒开始错动，边坡失稳破坏，此时黏结向右下方延伸破坏，颗粒间错动加剧。其中颗粒间黏结破坏形式以张拉破坏为主，当${F_{{\rm{FDM}}}}$=1.235时，张拉破坏近乎剪切破坏的2倍，边坡处于临界破坏状态。上述结果说明随着强度参数的不断劣化，由于边坡土体的抗拉强度弱，其破坏形式以张拉破坏为主，潜在滑动面贯通，边坡失稳破坏，位于滑动面处的土体产生较大位移，滑动面向边坡内部延伸发展。${F_{{\rm{FDM}}}}$=1.235为颗粒黏结破坏数量的突变点，可作为边坡的安全系数，与前文以特征点位移突变为判据耦合计算的安全系数相比，差距仅为2.59%，故现提出一种新的边坡失稳判据：颗粒黏结破坏数量突变判据。

图  11  不同强度折减系数下颗粒黏结状态云图

Figure  11.  Particle bonding state with different strength reduction factors

下载: 全尺寸图片 幻灯片

图  12  不同强度折减系数下颗粒黏结状态的演化过程

Figure  12.  Evolution of particle bonding state with different strength reduction factors

下载: 全尺寸图片 幻灯片

4   讨论

离散域位置、大小和形状的不同选择将对耦合强度折减法的计算精度和效率产生一定影响。因此，有必要讨论连续-离散耦合边坡模型中最优离散域的选取方法，为耦合强度折减法的精确和高效应用提供一定参考。总体上，离散域的选取应遵循以下的原则：

1)充分发挥离散元优势的原则。为了充分发挥离散元优势，应将离散域布置在应力和应变状态变化复杂的位置，比如边坡潜在滑动面区域，进而更有效地利用离散元分析边坡的细观力学变化规律。

2)权衡数值计算效率的原则。显然，离散域大小越大，单元数量越多，越能发挥离散元的优势，但这同时势必会降低数值计算的效率。已有研究表明，数值计算的时间与离散域内单元数量呈二次多项式增长关系^[30]。因此，确定离散域的大小应充分考虑数值模型的大小和计算机的运算能力，权衡数值计算效率。

3)简化数值计算模型的原则。离散域形状原则上应尽量规则、简单，一方面有助于连续-离散耦合墙区域的力和速度等信息交换的连续性和协调性，降低因形状畸变引起的模拟结果误差；另一方面亦便于研究人员的易操作性。

为了进一步提出最优离散域的选取方法，基于上述离散域的选取原则，设计离散域位置、大小和形状3个试验因素的数值模拟方案，开展耦合强度折减法计算精度研究，如表4所示。表中离散域相对位置是指离散域左上角顶点到坡面最短距离x与O点到坡面距离L的比值，如图13所示，规定x在坡面左侧为负，在坡面右侧为正；L可由式(3)计算；离散域相对大小是指离散域面积S与图13中三角形面积$S' = 0.5hl$的比值；离散域长宽比是指离散域长度a与离散域宽度b的比值：

表  4  离散域位置、大小和形状影响的数值模拟方案

Table  4.  Numerical simulation scheme for the influence of position, size and shape of discrete domain

试验分组		改变的试验因素	相对位置x/L	相对大小S/S′	相对形状a/b
大组	小组
Ⅰ	1	相对位置	−0.14	0.50	1.00
	2		0.20	0.50	1.00
	3		0.30	0.50	1.00
	4		0.35	0.50	1.00
	5		0.43	0.50	1.00
	6		0.46	0.50	1.00
	7		0.50	0.50	1.00
	8		0.65	0.50	1.00
Ⅱ	9	相对大小	0.43	0.18	1.00
	10		0.43	0.32	1.00
	11		0.43	0.50	1.00
	12		0.43	0.72	1.00
	13		0.43	0.98	1.00
Ⅲ	14	相对形状	0.43	0.32	0.50
	15		0.43	0.32	0.67
	16		0.43	0.32	1.00
	17		0.43	0.32	1.50
	18		0.43	0.32	2.00

下载: 导出CSV 

| 显示表格

式中，h和l分别为坡高和坡长。

图  13  离散域位置示意图

Figure  13.  Diagram of discrete region location

下载: 全尺寸图片 幻灯片

4.1   离散域位置引起的误差分析

图14为不同离散域位置下安全系数的变化。由图14可知，离散域相对位置x/L较小时，离散域边界过于靠近边坡坡面，因其边界效应对坡面附近区域的应力分布造成影响，离散域内颗粒黏结破坏范围明显增加，边坡安全系数波动较大，此时离散域位置对安全系数影响较大，且对边坡潜在滑动面范围造成一定的影响。随着离散域边界离坡面的距离增加，安全系数逐渐趋于标准值${F_{{\rm{FDM}}}}$=1.227，潜在滑动面贯通。其中F点和G点位置处离散域的滑动面未完全贯通，这是由于“切线法”计算安全系数所带来的误差，此时边坡仍处于临界破坏状态。以有限差分强度折减法计算出的安全系数${F_{{\rm{FDM}}}}$=1.227为标准，综合考虑安全系数的影响和潜在滑动面是否贯通，得到离散域最优位置x/L=0.43，其对应的安全系数为1.2229，差距仅为0.33%，此时边坡潜在滑动面已完全贯通。综上所述，选取离散域位置时，应保证边坡潜在滑动面穿过离散域，同时不宜离坡面太近。

图  14  不同离散域位置下安全系数的变化

Figure  14.  Variation of safety factor with different discrete region positions

下载: 全尺寸图片 幻灯片

4.2   离散域大小引起的误差分析

为消除离散域位置的影响，选取上述离散域最优位置，固定离散域左上角顶点位置，依照方案设置离散域的大小。图15为不同离散域大小下安全系数的变化。由图15可知，边坡安全系数变化较小并趋于标准值，边坡潜在滑动面均穿过离散域且贯通。其中当离散域相对大小S/S'=0.18时，离散域宽度小于边坡潜在滑动面竖向宽度，离散域范围偏小，且耦合边界处出现了一定数量的破坏。这说明在消除离散域位置的影响后，离散域大小对安全系数影响较小，但当离散域过小时，离散域不能完全覆盖边坡潜在滑动面所在范围，又由于受到尺寸效应的影响，其耦合边界处会出现一定程度的误差。当离散域过大时，数值计算的时间会随着颗粒数量的增加而成倍增加^[30-31]，极大地降低了计算效率。综合考虑边坡潜在滑动面范围、尺寸效应和计算效率的影响，得到离散域最优大小S/S'=0.32，其对应的安全系数为1.2229，差距仅为0.33%。综上所述，选取离散域大小时，应在充分发挥离散元优势的前提下，尽可能减少离散域的大小，有效地提高计算效率。

图  15  不同离散域大小下安全系数的变化

Figure  15.  Variation of safety factor with different discrete region sizes

下载: 全尺寸图片 幻灯片

4.3   离散域形状引起的误差分析

为消除离散域位置和大小的影响，选取上述离散域最优位置和大小，固定离散域左上角顶点位置，依照方案设置离散域的形状。图16为不同离散域形状下安全系数的变化。由图16可知，边坡安全系数变化较小并趋于标准值，边坡潜在滑动面均穿过离散域且贯通。离散域形状主要影响边坡潜在滑动面与离散域的相对位置关系。其中当离散域长宽比a/b=0.50时，离散域长度刚好为边坡潜在滑动面水平宽度；当a/b=1.50时，离散域宽度小于边坡潜在滑动面竖向宽度，离散域范围偏小。这说明在消除离散域位置和大小的影响后，离散域形状对安全系数影响较小，且对边坡潜在滑动面范围影响小。综上所述，综合考虑边坡潜在滑动面与离散域相对位置的影响，得到离散域最优长宽比a/b=0.67(1∶1.5)，其对应的安全系数为1.2267，差距仅为0.02%。综上所述，选取离散域形状时，应充分考虑边坡潜在滑动面与离散域的相对位置关系，尽量选择规则和简单的矩形，然后通过调整其长宽比来适应各种实际边坡工程。

图  16  不同离散域形状下安全系数的变化

Figure  16.  Variation of safety factor with different discrete region shapes

下载: 全尺寸图片 幻灯片

5   结论

本文开发出了一种连续-离散耦合强度折减法，基于某经典边坡算例验证了方法的可行性，并从细观角度上研究了边坡破坏机理和失稳判据，最后讨论了耦合强度折减法最优离散域的选取原则和方法。其主要结论如下：

(1)提出了宏、细观强度折减系数关系的建立方法，基于某经典算例计算发现宏观强度折减系数与细观强度折减系数并非呈线性相关，而是满足指数增长关系。

(2)耦合边界位移连续性很好，出现明显的分层现象，耦合强度折减法计算的安全系数与其他方法相比差距最大为5.39%，边坡潜在滑动面完全贯通，其位置、形状与其他方法相似度高，从多个角度验证了耦合强度折减法的可行性。

(3)边坡发生失稳破坏时，其破坏形式主要为张拉破坏，拉力链数量和颗粒间黏结破坏数量会出现明显的突变，由此提出2种新的边坡细观失稳判据：拉力链数量突变判据和颗粒黏结破坏数量突变判据。

(4)当离散域位置为x/L=0.425、大小为S/S′=0.32、形状为a/b=0.67时，耦合效果和计算效率最佳，其安全系数较标准值差距仅为0.02%，即为最优离散域。

(5)选取离散域时，应遵循充分发挥离散元优势、权衡数值计算效率、简化数值计算模型等原则，离散域应尽可能穿过边坡潜在滑动面，尽可能减少离散域的大小以提高计算效率，离散域形状宜采用规则和简单的矩形。