参考了陈育民老师的书，稍微修改了一点代码，发现在5.0里面运行后一直报错，每次运行到设置饱和密度的时候就会提示“system exception，forcing shutdown”
下面是代码
n
;===================================
;导入网格数据并生成“辅助”界面单元
impgrid water
interface 1 face ran x -0.1 100.1 y -0.1 100.1 z -0.1 100.1
set grav 0 -9.81 0
water den 1000
;===================================
;定义生成潜水面的函数
def water_table
p_i=i_head
p_ie=i_elem_head(p_i)
loop while p_ie # null
;返回三个相邻界面单元的三个顶点的地址
p_gp1=ie_vert(p_ie,1)
p_gp2=ie_vert(p_ie,2)
p_gp3=ie_vert(p_ie,3)
;将顶点坐标赋予给网格节点
x1=in_pos(p_gp1,1)
y1=in_pos(p_gp1,2)
z1=in_pos(p_gp1,3)
x2=in_pos(p_gp2,1)
y2=in_pos(p_gp2,2)
z2=in_pos(p_gp2,3)
x3=in_pos(p_gp3,1)
y3=in_pos(p_gp3,2)
z3=in_pos(p_gp3,3)
;以这三个节点，生成潜水面
command
water table face @x1,@y1,@z1 @x2,@y2,@z2 @x3,@y3,@z3
endcommand
p_ie=ie_next(p_ie)
endloop
end
;===================================
;进行边坡分析
impgrid kuaiti
group soil range group 1
group rock range group 2
;生成潜水面
@water_table
;删除为“辅助”网格单元和界面单元
interface 1 dele
dele zone range group bottom
dele zone range group top
;初始化材料密度: 饱和的和非饱和的
;施加边界约束条件
fix x y z rang z -0.1 0.1
fix x range x 99.9 100
fix x range x -0.1 0.1
fix y range y -0.1 0.1
fix y range y 99.9 100
;通过弹性求解生成初始应力场
model elas
def ini_dens
p_z = zone_head
loop while p_z # null
if z_group(p_z) = 'soil' then
if z_pp(p_z) # 0.0 then
z_density( p_z) = 2120
else
z_density( p_z) = 2020
endif
endif
if z_group( p_z) = 'rock' then
if z_pp( p_z) # 0.0 then
z_density(p_z) = 3460
else
z_density(p_z) = 3360
endif
endif
p_z = z_next(p_z)
endloop
end
@ini_dens
pro bulk 1e10 she 1e10 range group soil
pro bulk 1e10 she 1e10 range group rock
;位移和速度归零
ini xdisp 0 ydisp 0 zdisp 0
ini xvel 0 yvel 0 zvel 0
;定义材料特性
;冰碛土物理力学参数指标
model mohr
pro bulk 3.0e8 she 1.8e8 co 2.10e5 fric 25 ten 1e4 dil 10 range group soil
;强节理化辉长岩物理力学参数指标
pro bulk 8e8 she 6e8 co 4.60e5 fric 37 ten 4e3 dil 15 range group rock
hist n=5
hist unbal
plo hist 1
plo con zdis
step 6000
