求一个FLAC3D点安全系数命令流

边坡安全系数？flac3d里直接用solve fos就能算

这是Drucker Prager屈服条件的
new
rest 123.f3sav ;计算所得的最后一步结果
dele zone ran mo null ;删除开挖单元
def safety_factor ;定义安全系数
p_z=zone_head ;zone_head 指向单元体表
i=0
array z_coef(10000) ;定义数组
loop while p_z # null ;循环
i=i+1
o_coh=z_prop(p_z,'cohesion') ;得到计算中的 c 值
o_fri=z_prop(p_z,'friction') ;得到计算中的 fai 值
p_ten=z_prop(p_z,'tension') ;抗拉强度值
o_major=-1*z_sig1(p_z) ;计算中 sig1 值
o_middel=-1*z_sig2(p_z) ;计算中的 sig2 值
o_minor=-1*z_sig3(p_z) ;计算中的 sig3 值
o_ten=z_sig3(p_z)
II=o_major+o_middel+o_minor
JJ=((o_major-o_middel)^2+(o_middel-o_minor)^2+(o_minor-o_major)^2)/6
cosf=cos(o_fri*pi/180)
sinf=sin(o_fri*pi/180)
alf=(2*sinf)/(sqrt(3)*(3-sinf))
K=(6*o_coh*cosf)/(sqrt(3)*(3-sinf))
f_coef=abs((K-alf*II)/(sqrt(JJ))) ;得到安全系数
zone_id=z_id(p_z)
if o_ten>p_ten then
f_coef=0 ;材料拉破坏，用 0 做标识
endif
o_xpos=z_xcen(p_z)
o_ypos=z_ycen(p_z)
o_zpos=z_zcen(p_z)
z_coef(i)=string(o_xpos)+' '+string(o_ypos)+' '+string(o_zpos)+' '+string(f_coef)
p_z=z_next(p_z)
endloop
status=open('safety_factor.txt',1,1)
status=write(z_coef,i)
status=close
end
@safety_factor