这个程序用来模拟盾构隧道的开挖，得到地面的最终沉降值太小，不知道是什么原因，请各位大侠看看出了什么问题
开挖的过程如下：
new
gen zone cyl p1 3 0 0 p2 0 36 0 p3 0 0 3 size 3 30 6
gen zone cyl p1 0 0 -3 p2 0 36 0 p3 3 0 0 size 3 30 6
group tunnel
gen zone radcyl size 4 30 6 2 p1 4 0 0 p2 0 36 0 p3 0 0 4 dim 3 3 3 ratio 1 1 1 1.2
gen zone radcyl size 4 30 6 2 p1 0 0 -4 p2 0 36 0 p3 4 0 0 dim 3 3 3 ratio 1 1 1 1.2
gen zone brick p0 0 0 4 p1 add 4 0 0 p2 add 0 36 0 p3 add 0 0 14 size 3 30 8 ratio 1 1 1.1
gen zone brick p0 0 0 -4 p1 add 0 0 -9 p2 add 0 36 0 p3 add 4 0 0 size 6 30 3 ratio 1.1 1 1
gen zone brick p0 4 0 4 p1 add 20 0 0 p2 add 0 36 0 p3 add 0 0 14 size 10 30 8 ratio 1.1 1 1.1
gen zone brick p0 4 0 -13 p1 24 0 -13 p2 4 36 -13 p3 4 0 -4 p4 24 36 -13 p8 4 36 -4 p6 24 0 -4 p7 24 36 -4 size 10 30 6 ratio 1.1 1 0.909090909091
gen zone brick p0 4 0 0 p1 add 20 0 0 p2 add 0 36 0 p3 add 0 0 4 size 10 30 3 ratio 1.1 1 1
gen zone brick p0 4 0 -4 p1 24 0 -4 p2 4 36 -4 p3 4 0 0 p4 24 36 -4 p8 4 36 0 p6 24 0 0 p7 24 36 0 size 10 30 3 ratio 1.1 1 1
group rock range group tunnel not
mo mohr
pro bulk 8.3e6 she 3.85e6 fric 20 coh 18e3 ten 0 dil 0 range z -13 18
fix x range x -.1 .1
fix x range x 23.9 24.1
fix z range z -13.1 -12.9
fix y range y -.1 .1
fix y range y 35.9 36.1
set grav 0 0 -10
ini density 1800
ini szz -324e3 grad 0 0 18000
ini sxx -220e3 grad 0 0 12240
ini syy -220e3 grad 0 0 12240
hist unbal
hist id 2 gp zdisp 0 15 18
hist id 3 gp zdisp 2 15 18
hist id 4 gp zdisp 4 15 18
hist id 5 gp zdisp 6 15 18
hist id 6 gp zdisp 8 15 18
hist id 7 gp zdisp 10 15 18
hist id 8 gp zdisp 12 15 18
hist id 9 gp zdisp 14 15 18
hist id 10 gp zdisp 16 15 18
hist id 11 gp zdisp 18 15 18
hist id 12 gp zdisp 20 15 18
hist id 13 gp zdisp 22 15 18
hist id 14 gp zdisp 24 15 18
def sup_stress
lumda=0.9
o_press=-220000
o_grad=12240
s_press=o_press*lumda
s_grad=o_grad*lumda
end
sup_stress
plot contour zdisp
plot add sel geom scale=0.005 shrink=0.3 black black
plot set rotation 0 0 60
def excav
y0=-1.2
y0=y0+1.2
cut_i=y0/1.2+1
loop cut(cut_i,31)
y0=1.2*(cut-1)
y1=y0+1.2
yp0=1.2*cut-0.01
yp1=1.2*cut+0.01
yc0=y0-1.21
yc1=y0+0.01
command
mod null range group tunnel y y0 y1
end_command
command
apply nstress s_press grad 0 0 s_grad range cyl end1 0 yp0 0 end2 0 yp1 0 rad 3.0
end_command
if cut>1 then
command
sel shell id y0 range group rock cyl end1 0 yc0 0 end2 0 yc1 0 rad 3.0
sel node local xdir 1 0 0 ydir 0 -1 0 range x -0.1 0.1
sel node fix lsys range x -0.1 0.1
sel node fix x yr zr range x -0.1 0.1
sel node local xdir 0 0 -1 ydir 0 -1 0 range y -0.1 0.1
sel node fix lsys range y -0.1 0.1
sel node fix y xr zr range y -0.1 0.1
sel shell id y0 prop iso=(2.61e9,0.25) thick=0.3 density 2400
end_command
end_if
command
solve
end_command
if cut=5 then
command
save m1_6.sav
end_command
end_if
if cut=31 then
command
save m1_36.sav
end_command
end_if
end_loop
end
excav
ret
