求解围岩特征曲线

import itasca as itimport matplotlib.pyplot as mplit.command("python-reset-state false")# rock porperties of grade IV in N-m-s unit systemdens4 = 2600.0E4 = 2.0e9v4 = 0.3coh4 = 0.15e6fric4 = 30.0# lining porperties in N-m-s unit systemdensC = 2500.0EC = 31.5e9vC= 0.2cohC = 3.2e6fricC = 50.0disp = []zstress = []for i in range(4):    print(i)    buried_depth_inc = i*50    overburden = -1*dens4*9.81*buried_depth_inc    print(overburden)    it.command("""       model new       model largestrain off       zone import 'mesh.f3grid'       ;赋予本构及参数       zone cmodel assign mohr-coulomb       zone property young {} poisson {} density {} cohesion {} friction {}       ;赋予边界条件       zone face skin       zone face apply velocity-normal 0 range group 'west' or 'east'       zone face apply velocity-normal 0 range group 'south' or 'north'       zone face apply velocity-normal 0 range group 'bottom'       ;设置重力并求解初始地应力       model gravity 9.81       zone initialize-stresses overburden {}       zone face apply stress-normal {}       model solve elastic       ;清除初始位移、速度及塑性区       zone gridpoint initialize displacement (0,0,0)       zone gridpoint initialize velocity (0,0,0)       zone initialize state 0       ;开挖及支护       zone cmodel assign null range  group 'cut' or 'lining'       zone cmodel assign mohr-coulomb range group 'lining'       zone property young {} poisson {} density {} cohesion {} friction {} range group 'lining'       ;求解       model solve   """.format(E4,v4,dens4,coh4,fric4,overburden,overburden,EC,vC,densC,cohC,fricC))    gpnt = it.gridpoint.near((25,0.5,59.7))    zpnt = it.zone.find(6686)    disp.append(-1.0*gpnt.disp_z())    zstress.append(-1.0*zpnt.stress()[2][2])# 绘制围岩特征曲线mpl.plot(disp,zstress,color = 'black' ,linewidth = 1.0,marker = '*',markerfacecolor = 'red')mpl.title("Characteristic curve of surrounding rock")mpl.ylabel("stress/(Pa)")mpl.xlabel("disp/(m)")#C:\Users\18190\Desktop\case-01\为项目路径，根据实际路径进行替换mpl.savefig(r'C:\Users\18190\Desktop\case-01\buried_depth.png', dpi=600)mpl.show()

import itasca as itimport matplotlib.pyplot as mplit.command("python-reset-state false")# rock porperties of grade IV in N-m-s unit systemdens4 = 2600.0E4 = 2e9v4 = 0.3coh4 = 0.15e6fric4 = 30.0# lining porperties in N-m-s unit systemdensC = 2500.0EC = [31.5e10,31.5e9,31.5e8,31.5e7,31.5e6]vC= 0.2cohC = 3.2e6fricC = 50.0disp = []zstress = []for i in range(5):    print(i)    print(EC[i])    it.command("""        model new        model largestrain off        zone import 'mesh.f3grid'        ;赋予本构及参数        zone cmodel assign mohr-coulomb        zone property young {} poisson {} density {} cohesion {} friction {}        ;赋予边界条件        zone face skin        zone face apply velocity-normal 0 range group 'west' or 'east'        zone face apply velocity-normal 0 range group 'south' or 'north'        zone face apply velocity-normal 0 range group 'bottom'        ;设置重力并求解初始地应力        model gravity 9.81        zone initialize-stresses        model solve elastic        ;清除初始位移、速度及塑性区        zone gridpoint initialize displacement (0,0,0)        zone gridpoint initialize velocity (0,0,0)        zone initialize state 0        ;开挖及支护        zone cmodel assign null range group 'lining' or 'cut'        zone cmodel assign mohr-coulomb range group 'lining'        zone property young {} poisson {} density {}  cohesion {} friction {} range group 'lining'        ;求解        model solve    """.format(E4,v4,dens4,coh4,fric4,EC[i],vC,densC,cohC,fricC))    gpnt = it.gridpoint.near((25,0.5,59.7))    zpnt = it.zone.find(6686)    disp.append(abs(gpnt.disp_z()))    zstress.append(-1.0*zpnt.stress()[2][2])# 绘制围岩特征曲线mpl.plot(disp,zstress,color = 'r' ,linewidth = 1.0,marker = '*',markerfacecolor = 'black')mpl.title("Characteristic curve of surrounding rock")mpl.xlabel("disp/(m)")mpl.ylabel("stress/(Pa)")#路径说明同上mpl.savefig(r'C:\Users\18190\Desktop\case-01\E_changed.png', dpi=600)mpl.show()