在岩土工程数值模拟中，强度折减法（SRM）是计算边坡安全系数（FoS）最经典的方法之一。对于 FLAC3D 用户来说，计算安全系数通常只需要一句简单的内置命令：model factor-of-safety。这句指令强大且方便，但它像一个“黑盒”。底层到底是怎么折减的？如果我想在折减过程中加入自定义的力学准则（比如抗拉强度的顶点截断），或者想精准提取每一步折减的中间状态，内置命令往往显得不够灵活。今天，我们将借助 FLAC3D 强大的 Python (itasca 模块) 接口，手写一个基于二分法的自定义强度折减脚本，并让它与 FLAC3D 的内置 FoS 命令来一场“正面对决”！
自编强度折减法的Python 脚本主要分为四个核心步骤。下面详细拆解其中的几大细节：
第一步：建模与初始平衡态获取
任何折减的前提都是一个纯净的初始应力场。在赋予 Mohr-Coulomb 模型和重力后，需要求解初始平衡。
第二步：Python 二分法与抗拉截断（核心）
我们设定了初始的安全系数上下限（下限 1.0 代表稳定，上限 3.0 代表破坏），利用 while 循环进行二分法逼近，收敛容差设定为 0.001。
在这个环节，程序做了一个非常关键的工程力学优化——抗拉强度同步折减与顶点截断：
t_max_apex = c_trial / math.tan(phi_trial_rad)
t_trial = min(t_trial, t_max_apex)
由于莫尔-库仑准则的限制，材料的抗拉强度不能超过屈服包络线在正应力轴上的交点（顶点）。代码中加入了这个判定，确保物理意义的绝对严密！
第三步：状态判定与循环
每次生成新的折减参数后，恢复 initial_state.sav，赋予新参数，开始计算。
通过 zone.mech.ratio() 提取最大不平衡力比率。如果 ratio 小于 1e-5，说明边坡稳定，将当前系数赋给下限，并保存稳定存档；否则说明边坡破坏，赋给上限，保存破坏存档。
第四步：内置 FoS 登场与结论输出
为了验证我们自编逻辑的准确性，脚本最后会再次恢复初始状态，调用 FLAC3D 内置的 model factor-of-safety 命令计算安全系数（保存为 neizhi.sav）。
最终，通过 FISH 与 Python 的交互，将两个结果同框打印，进行直观对比。
这里展示脚本中最精髓的 Python 二分法折减控制段落：
while (f_upper - f_lower) > tol:
f_mid = 0.5 * (f_lower + f_upper)
# 1. 计算试验参数
c_trial = c_orig / f_mid
phi_trial_rad = math.atan(math.tan(phi_orig_rad) / f_mid)
phi_trial_deg = math.degrees(phi_trial_rad)
# 2. 抗拉强度顶点截断优化
t_trial = t_orig / f_mid
if phi_trial_rad > 0:
t_max_apex = c_trial / math.tan(phi_trial_rad)
t_trial = min(t_trial, t_max_apex)
# 3. 恢复初始状态并试算
it.command(f"model restore 'initial_state.sav'")
it.command(f"zone property cohesion {c_trial} friction {phi_trial_deg} tension {t_trial}")
it.command(f"model solve mechanical cycles {max_steps} ratio 1e-5")
# 4. 获取比率并判定收敛
ratio = it.fish.get("current_ratio")
if ratio < 1e-5:
f_lower = f_mid # 稳定
else:
f_upper = f_mid # 破坏
运行整个脚本后，将会在控制台看到极其干脆利落的对比结果：
Custom Python Bisection FoS vs FLAC3D Built-in Command FoS。
通常情况下，两者结果会高度一致（误差在设定的 tol 范围内）。
下面是本人拿一个边坡模型进行的测试结果对比
老规矩，有需要的小伙伴们，欢迎关注、点赞、收藏、转发，完成后通过留言发送邮箱，只需要稍微修改收到的py文件，将模型替换自己的网格文件路径与材料参数，在 FLAC3D 中运行即可。
