1. 边坡安全系数的计算

边坡可以分为土质边坡和岩质边坡，如果是土质边坡，一般可以看成连续介质，使用有限元软件或者FLAC等进行计算，用3DEC/UDEC等离散元计算边坡的安全系数一般会考虑含节理的岩质边坡类型，当然也可以考虑不含节理的，下面分别介绍。

1.1 基于结构面接触的边坡安全系数计算

对于含有结构面的边坡模型，可以采用接触判断计算块体接触面的安全系数，计算原理是结构面的抗滑力与下滑力的比值作为安全系数的定义。

对于常用的摩尔库伦节理模型，结构面上的抗滑力主要由下面的公式计算:

下滑力用结构面上的切向力累加得到，具体实现可以参考下面的3DEC5.2的FISH函数，UDEC6.0可以进行类似的修改即可。

以下函数参考了《块体离散元数值模拟技术及工程应用》——石崇等编著一书

如果显示不全，点击代码框顶部可以全屏查看

def call_fos

            ci = c_near(50,-60,105)   ;获取结构面的接触索引，c_near中的坐标为该结构面上某一点
            ;或者使用接触的编号进行判断接触，前提是设置jset的时候定义了编号或者你知道该接触的编号
            ;_joint_ID = 10
            ;
            ;
            ;_ci = contact_head
            ;
            ;loop while _ci # 0
            ;
            ;_contact_id = c_id(_ci)
            ;
            ;if _contact_id == _joint_ID then

            _ci_sub = c_cx(ci)        ;获取该接触的子接触索引

            sumSubcontactNF = 0.0     ;用于累加子接触上的法向力
            sumSubcontactSF = 0.0     ;用于累加子接触上的切向力
            sumSubcontactNum = 0      ;用于累加子接触的数量
            sumSubcontactArea = 0.0   ;用于累加子接触的面积
            sumResistingForce=0.0     ;用于累加抗滑力
            sumDrivingForce=0.0       ;用于累加下滑力
            _fos=0.0                  ;用于存储安全系数

            loop while _ci_sub # 0                                        ;遍历该接触下的所有子接触

                subcontactNF = cx_nforce(_ci_sub)                         ;获取子接触的法向力
                sumSubcontactNF = sumSubcontactNF + subcontactNF          ;累加子接触上的法向力
                ;用于计算子接触的切向力，可以用分量计算也可以求矢量的模
                ;subcontactSF = sqrt((cx_xsforce(_ci_sub))^2+(cx_ysforce(_ci_sub))^2+(cx_zsforce(_ci_sub))^2)
                 subcontactSF = mag(cx_sforce(_ci_sub))
                sumSubcontactSF = sumSubcontactSF + subcontactSF          ;累加子接触上的法向力

                subcontactArea = cx_area(_ci_sub)                         ;获取子接触的面积
                sumSubcontactArea = sumSubcontactArea + subcontactArea    ;累加子接触上的面积

               sumSubcontactNum =  sumSubcontactNum + 1                   ;累加子接触上的数目

               _ci_sub = cx_next(_ci_sub)                                 ;获取下一子接触的索引

               jcoh_j_id = cx_prop(_ci_sub,'jcoh')                          ;获取结构面的粘聚力数值
               jfri_j_id = cx_prop(_ci_sub,'jfri')                          ;获取结构面的摩擦角数值
               ;计算并累加抗滑力
               sumResistingForce = sumResistingForce + jcoh_j_id * subcontactArea +subcontactNF  * tan(jfri_j_id * degrad)

            endloop                                                       ;结束循环

            sumDrivingForce = sumSubcontactSF                              ;切向力视为下滑力

            if sumDrivingForce # 0
            _fos = sumResistingForce/sumDrivingForce                       ;计算安全系数
            endif

            ;endif
            ;_ci =c_next(_ci)
            ;endloop
end

1.2 基于强度折减法的边坡安全系数计算

以下内容参考《FLAC/FLAC3D基础与工程实例第二版》——陈育民徐鼎平编著一书

3DEC和UDEC内置的强度折减法应该也是以二分法进行快速寻找求解的，以下命令流为FLAC3D 3.0代码实现的二分法求解强度折减系数的过程，对上一节所述的原理的理解有一定帮助，也对自编强度折减法有一定的启示作用。

new
;建立网格模型
gen zone brick & 
p0 0 0 0 pl 2 0 0 p2 0 0.5 0 p3 0 0 3 size 3 1 3 
gen zone brick & 
p0 2 0 0 pl 20 0 0 p2 2 0.5 0 p3 2 0 3 size 17 1 3 & 
ratio 1.03 1 1 
gen zone brick & 
pO 2 0 3 pl 20 0 3 p2 2 0.5 3 p3 12 0 13 & 
p4 20 0.5 3 p5 12 0.5 13 p6 20 0 13 & 
p7 20 0.5 13 size 17 1 17 ratio 1.03 1 1

;****************************************
;自定义二分法强度折减法
def SSR                                                 ;定义强度折减法函数
;定义有关参数及循环终止条件
aitl  = 0.02                                                    ;定义安全系数上下限之差的临界值（用户自定义）
k11   = 1.0                                                     ;初始安全系数下限值（用户自定义）
k12   = 2.0                                                     ;初始安全系数上限值（用户自定义）
ks    = (k11+k12)/2                                             ;折减系数的第一次更新

loop while (k12-k11)>aitl                                    ;定义计算循环的终止条件
    coh1  = 12380/ks                                            ;定义粘结力折减公式
    fri1  = (atan((tan(20*pi/180))/ks))* 180/pi                 ;定义摩擦角折减公式
    dila1 = 20.0                                                ;赋予剪胀角
    ten1  = 1e6                                                 ;赋予抗拉强度
    grav0 = -10                                                 ;赋予重力加速度
    dens1 = 2000                                                ;赋予密度
    K1    = 1e8                                                 ;赋予体积模量
    G1    = 3e7                                                 ;赋予剪切模量

;折减的实现过程                                                ;
    command                                                     ;开始执行下面的命令流

        model null                                              ;设置空模型

        ;初始应力场的生成                                       ;
        model elastic                                           ;设置弹性本构模型
        pro bulk 1e10 she 3e9 dens dens1                        ;设置物理力学参数
        fix x y z range z -0.1 0.1                              ;定模型底部边界的x、y、z 方向速度 
        fix x range x 19.9 20.1                                 ;固定模型边界x=20面所有点x 方向速度 
        fix x range x -0.1 0.1                                  ;固定模型边界x=0面所有点x 方向速度 
        fix y                                                   ;固定模型边界所有点y方向速度 
        set grav 0 0 grav0                                      ;设置重力加速度
        solve                                                   ;求解获得初始应力场
        ini xdisp 0 ydisp 0 zdisp 0                             ;位移场清零
        ini xvel 0 yvel 0 zvel 0                                ;速度场清零

        ;塑性阶段求解                                           ;设置摩尔-库仑本构模型
        model mohr pro bulk K1 she G1 dens dens1 coh coh1 &     ;设置物理力学参数
        friction fri1 dil dila1 tens ten1                       ;
        set mech ratio 9.8e-6                                   ;设置力不平衡比率
        solve step 10000                                        ;设置求解极限时步

    endcommand                                                  ;结束上一段命令的执行

    ;二分法的实现过程                                           ;
    if mech_ratio < 1.0e-5                                   ;如果当前折减系数下在规定步数内计算达到平衡
        k11 = ks                                                ;将当前折减系数赋给下限
        k12 = k12                                               ;上限值保持不变
    else                                                        ;如果前折减系数下在规定步数内计算未达到平衡
        k12 = ks                                                ;将当前折减系数赋给上限
        k11 = k11                                               ;下限值保持不变
    endif                                                       ;
                                                                ;
    ks= (k11+k12)/2                                             ;更新折减系数
                                                                ;
endloop                                                         ;重复循环计算

;计算结果的保存                                                ;定义结果保存的文件名
fosfile0 = '_fos'+'.sav'                                        ;  
command                                                         ;
    save fosfile0                                               
endcommand                                                      ;

end                                                             ;结束整个强度折减法函数的自定义过程

;****************************************                       ;
;程序执行及结果显示                                          ;
SSR                                                             ;执行自定义的强度折减法函数  
pr ks                                                           ;显示安全系数

下面的内容是对上述命令流的解释：

这里对这一段命令文件的思路做些说明，以便读者更容易理解。该命令文件中，model null 是关键，其目的是通过这一命令清除先前的计算结果，同时更新安全系数的上、下限值，然后采用其他非空本构模型来激活网格单元，接着进行下一次强度折减后的计算。如此反复循环计算，直至上、下限值之差小于0.02时，终止整个强度折减法的执行，此时上、下限的均值即为最终的安全系数。
1.力不平衡力比率R 小于临界值 9.8e10^-6 时；
2.运行时步超过10000时步时。
折减后的单次计算过程在满足上述任何一种情况时，便退出当前计算，开始新的折减计算过程。之所以不完全以R小于临界值作为单次折减计算过程的终止条件，是因为有时候会出现这种情况：体系可能会在力不平衡比率远小于9.8e10^-6时，便已达到稳定的塑性流动（实际上已达到力平衡），计算过程即使再进行下去，力不平衡比率也不可能达到这一临界值。对这种情况，该标准显然过于苛刻，也是不可能达到的，这时便可借助定义运行时步极限即第二条标准来实现计算过程的终止。
不过，运行时步极限往往需要事先粗略试算确定（即对应内置折减法中体系典型时步此的估算），该时步一般要大于体系达到力平衡状态时所对应时步的临界值。通过上述说明可以看出，这段自定义强度折减法的单一计算过程是以第一条标准为首选终止标准，只有当第一条终止标准无法达到时，才启用第二条标准，这实际上与FLAC3D内置强度折减法计算过程的终止标准是大体一致的，这有效确保了该自编强度折减法的通用性。