利用python实现二分法求解边坡安全系数

1年前浏览6128

摘要：利用python实现二分法求解边坡安全系数。

1.问题来源

本例为陈育民、徐鼎平主编《FLAC、FLAC3D基础与工程实例》（第二版）第十四章案例，案例详见下图：

2.代码

代码思路参考徐文刚,余旭荣,年廷凯等.基于FLAC3D的三维边坡稳定性强度折减法计算效率改进算法及其应用[J].吉林大学学报(地球科学版),2021,51(05):1347-1355.利用python自编强度折减法代码可避免重复计算自重应力场，有效缩短计算时长。

import itasca as itimport mathit.command("python-reset-state false")
G1 = 3e7K1 = 1e8dens1 = 2000.0grav = -10.0ten1 = 1e6dial1 = 20.0coh1 = 12380fri1 = 20.0it.command("""           model new           model largestrain off           zone create brick point 0 0 0 0 point 1 2 0 0 point 2 0 0.5 0 point 3 0 0 3 size 3 1 3           zone create brick point 0 2 0 0 point 1 20 0 0 point 2 2 0.5 0 point 3 2 0 3 size 17 1 3 ratio 1.03 1 1           zone create brick point 0 2 0 3 point 1 20 0 3 point 2 2 0.5 3 point 3 12 0 13 point 4 20 0.5 3 point 5 12 0.5 13 point 6 20 0 13 point 7 20 0.5 13 size 17 1 17 ratio 1.03 1 1           zone cmodel assign elastic           zone property bulk 1e10 shear 3e9 density {}            zone face skin break-angle 60           zone face apply velocity-normal 0 range group 'South' or 'North'           zone face apply velocity-normal 0 range group 'East' or 'West'           zone face apply velocity (0,0,0) range group 'Bottom'           model gravity 0,0,{}           model solve            zone gridpoint initialize displacement (0,0,0)           zone gridpoint initialize velocity (0,0,0)           model save 'ini.f3sav' """.format(dens1,grav))
ait1 = 0.02k11 = 1.0k12 = 2.0ks = (k11+k12)/2.0
while k12 - k11 > ait1:    coh1 = 12380.0/ks    fri1 = math.atan(math.tan(20.0*math.pi/180.0)/ks)*180.0/math.pi    dial1 = fri1    it.command("""               model restore 'ini.f3sav'               zone cmodel assign mohr-coulomb               zone property bulk {} shear {} density {} cohesion {} friction {} dilation {} tension {}               model solve ratio 9.8e-6 or cycles 10000""".format(K1,G1,dens1,coh1,fri1,dial1,ten1))    if it.zone.mech_ratio() < 1.0e-5:        k11 = ks        k12 = k12    else:        k12 = ks        k11 = k11    ks = (k11+k12)/2.0print(ks)