利用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.0while 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)