命令流
1 gen zone bri p0 0 0 0 p1 10 0 0 p2 0 10 0 p3 0 0 10 & p4 15 15 0 p5 0 15 15 p6 15 0 10 p7 20 20 20 & size 10 10 10 rat 1.0 0.9 1.1 group brick_1 gen zone bri p0 20 0 0 p1 add 10 0 0 p2 add 0 20 0 p3 add 0 0 15 & size 10 10 10 rat 1.0 0.9 1.1 group brick_2 gen zone bri p0 40 0 0 edge 10 size 10 10 10 rat 1.0 0.9 1.1 group brick_3 plot sur 2-1 new gen zon bri size 3 3 3 model elas prop bulk 3e8 shear 1e8 ini dens 2000 fix z ran z -.1 .1 fix x ran x -.1 .1 fix x ran x 2.9 3.1 fix y ran y -.1 .1 fix y ran y 2.9 3.1 set grav 0 0 -10 solve app nstress -10e4 ran z 3 x 1 2 y 1 2 hist gp vel 0 0 3 hist gp vel 0 3 3 plo hist 1 red plo add hist 2 blue solve 3-1 ;-------------------------------------工程信息
;Project Record Tree export ;Title:Simple test ;---------------------------------计算第一步
;... STATE: STATE1 .... config grid 10,10 model elastic group 'User:Soil' notnull model elastic notnull group 'User:Soil'
prop density=1500.0 bulk=3E6 shear=1E6 notnull group 'User:Soil' fix x y j 1 fix x i 1 fix x i 11 set gravity=9.81
history 999 unbalanced solve
save state1.sav ;----------------------------------计算第二步 ;... STATE: STATE2 .... initial xdisp 0 ydisp 0 initial xvel 0 yvel 0 model null i 4 7 j 8 10 group 'null' i 4 7 j 8 10 group delete 'null' history 1 xdisp i=4, j=11 solve save state2.sav ;--------------------------------绘图命令 ;*** plot commands **** ;plot name: syy plot hold grid syy fill ;plot name: Unbalanced force plot hold history 999 ;plot name: grid plot hold grid magnify 20.0 lred grid displacement ;plot name: Xdis-A plot hold history 1 line 5-1 new ; =============================== ; 定义球体半径和半径方向上单元网格数
; =============================== def parm rad=10.0
rad_size=5 end parm ; =============================== ; 建立八分之一球体外接立方体网格
; ===============================
gen zone pyramid p0 rad 0 0 p1 rad 0 rad p2 rad rad 0 p3 0 0 0 &
p4 rad rad rad size rad_size rad_size rad_size group 1 gen zone pyramid p0 0 rad 0 p1 rad rad 0 p2 0 rad rad p3 0 0 0 &
p4 rad rad rad size rad_size rad_size rad_size group 2 gen zone pyramid p0 0 0 rad p1 0 rad rad p2 rad 0 rad p3 0 0 0 &
p4 rad rad rad size rad_size rad_size rad_size group 3 ; ================================== ; 利用FISH语言将内部立方体节点调整到球面 ; ================================== def make_sphere p_gp=gp_head
loop while p_gp#null
; 获取节点点坐标值：P=(px,py,pz) px=gp_xpos(p_gp) py=gp_ypos(p_gp) pz=gp_zpos(p_gp)
dist=sqrt(px*px+py*py+pz*pz) if dist>0 then ; 节点位置调整
maxp=max(px,max(py,pz)) k=(maxp/rad)*(rad/dist) gp_xpos(p_gp)=k*px gp_ypos(p_gp)=k*py gp_zpos(p_gp)=k*pz end_if
p_gp=gp_next(p_gp) end_loop end
make_sphere
; =============================== ; 利用镜像生成完整球体网格
; =============================== gen zone ref gen zone ref dip 90 gen zone ref dip 90 dd 90
; =============================== ; 显示球体网格
; =============================== plot surf
pl set back wh pl bl gr 5-3 n
gen zon bri size 1 1 2 group soil ran z 1 1 group rock ran z 0 1 expgrid 1.flac3d 6-1
gen zon bri size 3 3 3 model mohr
prop bu 3e6 sh 1e6 coh 10e3 fric 15 fix z ran z -.1 .1 fix x ran x -.1 .1 fix x ran x 2.9 3.1 fix y ran y -.1 .1 fix y ran y 2.9 3.1 ini dens 2000 hist unbal set grav 10 solve elastic save 6-1.sav 6-2
rest 6-1.sav
ini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0
app nstress -100e3 ran z 2.9 3.1 x 1 2 y 1 2 solve
save 6-2.sav 6-3
rest 6-1.sav
ini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0
app nstress -100e3 ran z 2.9 3.1 x 1 2 y 1 2 hist id=2 gp zdis 1 1 3 hist id=3 gp zdis 1 1 2 hist id=4 gp xdis 1 1 3 hist id=5 gp xdis 1 1 3 hist id=6 zone szz 1 1 3
hist id=7 zone szz 1.5 1.5 2.5 hist id=8 zone sxz 1.5 1.5 2.5
solve
save 6-3.sav 6-4
rest 6-3.sav set log on
set logfile 6-2.log print zone stress print gp dis set log off 6-5
rest 6-1.sav
ini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0
app nstress -100e3 ran z 2.9 3.1 x 1 2 y 1 2 plot set rot 20 0 30
plot con szz ou on magf 10 plot add hist 1
set movie avi step 1 file 6-5.avi movie start solve
movie finish 7-1 n
gen zon bri size 1 1 2 model elas
prop bulk 3e7 shear 1e7 fix z ran z 0 fix x ran x 0 fix x ran x 1 fix y ran y 0 fix y ran y 1 ini dens 2000 set grav 0 0 -10 solve plo con sz 7-2 n
gen zon bri size 1 1 2 model mohr
prop bulk 3e7 shear 1e7 c 1e10 f 15 tension 1e10 fix z ran z 0 fix x ran x 0 fix x ran x 1 fix y ran y 0 fix y ran y 1 ini dens 2000 set grav 0 0 -10 solve
prop bulk 3e7 shear 1e7 c 10e3 f 15 ten 0 solve plo con sz 7-3 n
gen zone brick size 1 1 2 model mohr
prop bulk 3e7 shear 1e7 coh 10e3 fri 15 ten 0 fix z ran z 0 fix x ran x 0 fix x ran x 1 fix y ran y 0 fix y ran y 1 ini dens 2000 set grav 0 0 -10 solve elas plo con sz 7-4 new
gen zone brick size 1 1 2 model mohr
prop bulk 3e7 shear 1e7 coh 10e3 fri 15 ten 0 fix z ran z 0 fix x ran x 0 fix x ran x 1 fix y ran y 0 fix y ran y 1 ini dens 2000
ini szz -40e3 grad 0 0 20e3 ran z 0 2 ini syy -20e3 grad 0 0 10e3 ran z 0 2 ini sxx -20e3 grad 0 0 10e3 ran z 0 2
set grav 0 0 -10 solve plo con sz 7-5 n
gen zon bri size 1 1 2 model m
prop bulk 3e7 shear 1e7 c 10e10 f 15 ten 1e10 fix z ran z 0 fix x ran x 0 fix x ran x 1 fix y ran y 0 fix y ran y 1
ini dens 2000 ran z 0 1 ini dens 1500 ran z 1 2
ini szz -35e3 grad 0 0 20e3 ran z 0 1 ini syy -17.5e3 grad 0 0 10e3 ran z 0 1 ini sxx -17.5e3 grad 0 0 10e3 ran z 0 1 ini szz -15e3 grad 0 0 15e3 ran z 1 2 ini syy -7.5e3 grad 0 0 7.5e3 ran z 1 2 ini sxx -7.5e3 grad 0 0 7.5e3 ran z 1 2 ini pp 10e3 grad 0 0 -10e3 ran z 0 1 set grav 0 0 -10 solve plo con sz 7-6 n
gen zon bri size 1 1 2 model m
prop bulk 3e7 shear 1e7 c 10e10 f 15 ten 1e10 fix z ran z 0 fix x ran x 0 fix x ran x 1 fix y ran y 0 fix y ran y 1
ini dens 2000 ran z 0 2
ini szz -50e3 grad 0 0 20e3 ran z 0 1 ini syy -30e3 grad 0 0 10e3 ran z 0 1 ini sxx -30e3 grad 0 0 10e3 ran z 0 1 ini pp 30e3 grad 0 0 -10e3 ran z 0 2 app nstress -10e3 ran z 2 set grav 0 0 -10 solve plo con sz 7-7 new
gen zone brick p0 0 0 0 p1 60 0 0 p2 0 60 0 p3 0 0 90 & p4 60 60 0 p5 0 60 90 p6 60 0 150 p7 60 60 150 & size 6 6 10 model elas
pro bulk 10e10 she 10e10 ini den 2500
apply sxx -1e9 grad 0 0 1.1111111e7 range x -.1 .1 apply sxx -1e9 grad 0 0 6.6666666e6 range x 59.9 60.1 apply syy -1e9 grad 0 0 8.3333333e6 range y -.1 .1 apply syy -1e9 grad 0 0 8.3333333e6 range y 59.9 60.1 apply szz -1e8 grad 0 0 8.3333333e5 ran z 0 120 set grav 0 0 -10 step 30000
ini xdisp 0 ydisp 0 zdisp 0 ini xvel 0 yvel 0 zvel 0 plo cont szz 7-8 new
gen zone brick p0 0 0 0 p1 60 0 0 p2 0 60 0 p3 0 0 90 & p4 60 60 0 p5 0 60 90 p6 60 0 150 p7 60 60 150 & size 6 6 10 model elas
pro bulk 10e10 she 10e10 ini den 2500
ini sxx -1e9 grad 0 0 1.1111111e7 range x -.1 .1 ini sxx -1e9 grad 0 0 6.6666666e6 range x 59.9 60.1 ini syy -1e9 grad 0 0 8.3333333e6 range y -.1 .1 ini syy -1e9 grad 0 0 8.3333333e6 range y 59.9 60.1 ini szz -1e8 ran z -.1 .1 fix x y z ran z -.1 .1 set grav 0 0 -10 solve
ini xdisp 0 ydisp 0 zdisp 0 ini xvel 0 yvel 0 zvel 0 plo cont szz
7-9 new
gen zone brick p0 0 0 -50 p1 27.5 0 -50 p2 0 5 -50 p3 0 0 -10 size 8 1 10 group clay
gen zone brick p0 27.5 0 -50 p1 100 0 -50 p2 27.5 5 -50 p3 27.5 0 -10 ratio 1.1 1 1 size 12 1 10 group clay
gen zone brick p0 0 0 -10 p1 27.5 0 -10 p2 0 5 -10 p3 0 0 0 ratio 1 1 0.8 size 8 1 4 group soil
gen zone brick p0 27.5 0 -10 p1 100 0 -10 p2 27.5 5 -10 p3 27.5 0 0 ratio 1.1 1 0.8 size 12 1 4 group soil
gen zone brick p0 0 0 0 p1 27.5 0 0 p2 0 5 0 p3 0 0 5 p4 27.5 5 0 &
p5 0 5 5 p6 20 0 5 p7 20 5 5 size 8 1 5 group dam
fix x y z ran z -49.9 -50.1 fix x ran x -.1 .1
fix x ran x 99.9 100.1 fix y
model mohr ran z -50 0 model null ran z 0 5
prop bulk 7.8e6 shear 3.0e6 coh 10e10 tension 1e10 ran group soil
ini dens 1500 ran group soil
prop bulk 3.91e6 shear 1.5e6 coh 10e10 tension 1e10 ran group clay
ini dens 1800 ran group clay
set grav 0 0 -9.8 hist id=1 unbal solve
prop bulk 7.8e6 shear 3.0e6 coh 10e3 fric 15 ran group soil
prop bulk 3.91e6 shear 1.5e6 coh 20e3 fric 20 ran group clay solve
save elastic.sav
ini xdis 0 ydis 0 zdis 0 ;将节点位移清零 ini xvel 0 yvel 0 zvel 0 ;将节点速度清零
hist id=2 gp zdis 0 0 0 ;记录地基顶部中心点的沉降 hist id=3 gp zdis 27.5 0 0 ;记录路基坡脚处的沉降 hist id=4 gp xdis 27.5 0 0 ;记录路基坡脚处的水平
位移
model elastic ran z 0 1 ; ;激活0 m ~ 1 m的单元 prop bulk 7.8e6 shear 3.0e6 ran z 0 1 ini dens 1500 ran z 0 1 solve ;按软件默认精度求解 save fill-1.sav
model elastic ran z 1 2
prop bulk 7.8e6 shear 3.0e6 ran z 1 2 ini dens 1500 ran z 1 2 solve
save fill-2.sav
model elastic ran z 2 3
prop bulk 7.8e6 shear 3.0e6 ran z 2 3 ini dens 1500 ran z 2 3 solve
save fill-3.sav
model elastic ran z 3 4
prop bulk 7.8e6 shear 3.0e6 ran z 3 4 ini dens 1500 ran z 3 4 solve
save fill-4.sav
model elastic ran z 4 5
prop bulk 7.8e6 shear 3.0e6 ran z 4 5 ini dens 1500 ran z 4 5 solve
save fill-5.sav pau
;plo bl gr
;gen zone brick p0 0 0 0 p1 100 0 0 p2 0 5 0 p3 0 0 5 size
gen zone brick p0 0 0 -50 p1 27.5 0 -50 p2 0 5 -50 p3 0 0 -10 size 8 1 10 group clay
gen zone brick p0 27.5 0 -50 p1 100 0 -50 p2 27.5 5 -50 p3 27.5 0 -10 ratio 1.1 1 1 size 12 1 10 group clay
gen zone brick p0 0 0 -10 p1 27.5 0 -10 p2 0 5 -10 p3 0 0 0 ratio 1 1 0.8 size 8 1 4 group soil
gen zone brick p0 27.5 0 -10 p1 100 0 -10 p2 27.5 5 -10 p3 27.5 0 0 ratio 1.1 1 0.8 size 12 1 4 group soil gen zone brick p0 0 0 0 p1 27.5 0 0 p2 0 5 0 p3 0 0 5 p4 27.5 5 0 &
p5 0 5 5 p6 20 0 5 p7 20 5 5 size 8 1 5 group dam 7-10
set log on ;打开log记录
set logfile 1.log ;设置记录文件名为：1.log restore fill-1.sav ;调用保存的文件 print gp dis range id 517 any id 533 any ;输出两个节点的变形值
restore fill-2.sav
print gp dis range id 517 any id 533 any restore fill-3.sav
print gp dis range id 517 any id 533 any restore fill-4.sav
print gp dis range id 517 any id 533 any restore fill-5.sav
print gp dis range id 517 any id 533 any set log off ;关闭log记录 8-1 def abc
abc = 1 + 2 * 3 abcd = 1.0 / 2.0 end abc
print fish 8-2 def abc
if aa < 0 then abc = 0.0 else
abc = 2.0 * aa endif end abc 8-3
def abc
loop aa (1, 2.5) command print aa endcommand endloop end 8-4 new
gen zon bri size 3 3 3 model elastic
prop bu 3e7 sh 1e7 ini dens 2000
fix x y z ran z -.1 .1 fix x ran x -.1 .1 fix x ran x 2.9 3.1 fix y ran y -.1 .1 fix y ran y 2.9 3.1 set grav 10 solve
ini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0 save 8-4.sav 8-5
rest 8-4.sav def E_modify
p_z = zone_head d_k = 704 d_n = 0.38
d_pa = 101325.0 ;//标准大气压 loop while p_z # null
sigma_3 = -1.0 * z_sig1(p_z)
E_new = d_k * d_pa * (sigma_3 / d_pa) ^ d_n z_prop(p_z,'young') = E_new p_z = z_next(p_z) endloop end
E_modify 8-6
rest 8-5.sav
table 1 name load_settlement def add_load
p_gp = gp_near(2,1,3) loop n (1,5)
app_load = n * (-1000e3)
file_name = '7-6_add_step' + string(n) + '.sav' command
app nstress app_load ran z 2.9 3.1 x 1 2 y 1 2 solve
save file_name endcommand
xtable(1,n) = -1.0 * app_load ytable(1,n) = gp_zdisp(p_gp) endloop end
add_load save 8-6.sav 8-7
rest 8-6.sav
def find_max_disp p_gp = gp_head maxdisp_value = 0.0 maxdisp_gpid = 0 loop while p_gp # null
disp_gp = sqrt(gp_xdisp(p_gp) ^ 2 + gp_ydisp(p_gp) ^ 2 + gp_zdisp(p_gp) ^ 2)
if disp_gp > maxdisp_value maxdisp_value = disp_gp
maxdisp_gpid = gp_id(p_gp) endif
p_gp = gp_next(p_gp) endloop end
find_max_disp
print maxdisp_value maxdisp_gpid
rest 8-6.sav config zextra 1 def get_sigma_dif p_z = zone_head loop while p_z # null
sigma_dif = z_sig3(p_z) - z_sig1(p_z) z_extra(p_z,1) = sigma_dif p_z = z_next(p_z) endloop end
get_sigma_dif plot con zextra 1 9-1 ;--------------------------------------------------- ;
; 移来移去法接触面的建立
;--------------------------------------------------- n
gen zone radcyl p0 (0,0,0) p1(8,0,0) p2 (0,0,-5) p3 (0,8,0) &
p4 (8,0,-5) p5 (0,8,-5) p6 (8,8,0) p7 (8,8,-5) & p8 (.3,0,0) p9 (0,.3,0) p10 (.3,0,-5) p11 (0,.3,-5) & size 3 10 6 15 ratio 1 1 1 1.15
gen zone radcyl p0 (0,0,-5) p1 (8,0,-5) p2 (0,0,-8) p3 (0,8,-5) &
p4(8,0,-8) p5 (0,8,-8) p6 (8,8,-5) p7 (8,8,-8) &
p8 (.3,0,-5) p9 (0,.3,-5) p10 (.3,0,-8) p11 (0,.3,-8) & size 3 6 6 15 ratio 1 1 1 1.15 fill gen zone reflect dd 270 dip 90 group clay ;
interface 1 face range cylinder end1 (0,0,0) end2 (0,0,-5.1) radius .31 &
cylinder end1 (0,0,0) end2 (0,0,-5.1) radius .29 not interface 2 face range cylinder end1 (0,0,-4.9) end2 (0,0,-5.1) radius .31 ; pause
gen zone cyl p0 (0,0,6) p1 (.3,0,6) p2 (0,0,1) p3 (0,.3,6) &
p4 (.3,0,1) p5 (0,.3,1) & size 3 10 6
gen zone cyl p0 (0,0,6.1) p1 (.3,0,6.1) p2 (0,0,6) p3 (0,.3,6.1) &
p4 (.3,0,6) p5 (0,.3,6) & size 3 1 6
gen zone reflect dd 270 dip 90 range z 1 6.1 group pile range z 1 6.1 pause
ini z add -6.0 range group pile save pile_geom.sav 9-2
;导来导去法
;--------------------------------------------- n
gen zone brick size 3 3 3
group 2 range x 1 2 y 1 2 z 1 2 group 1 range gr 2 not save 1.sav
del ran group 2 not
interface 1 face range x 1 y 1 2 z 1 2 interface 1 face range x 2 y 1 2 z 1 2 interface 1 face range x 1 2 y 1 z 1 2 interface 1 face range x 1 2 y 1 2 z 1 interface 1 face range x 1 2 y 1 2 z 2
rest 1.sav
del ran group 2 expgrid 1.fac3d
impgrid 1.flac3d
model ela
fix x y z ran z 0 ini den 2000 set grav 0 0 -10
interface 1 prop kn 20e6 ks 20e6 coh 10e3 fri 15 app nstr -200e3 ran x 0 1 y 1 2 z 3 solve 9-3
;切割模型法
;------------------------------ n
gen zone brick size 3 3 3
group 1 range x 1 2 y 1 2 z 2 3 group 2 range group 1 not gen separate 1 int 1 wrap 1 2 int 1 maxedge 0.5 plo int red 9-4
; Create Material Zones gen zone brick size 5 5 5 &
p0 (0,0,0) p1 (3,0,0) p2 (0,3,0) p3 (0,0,5) & p4 (3,3,0) p5 (0,5,5) p6 (5,0,5) p7 (5,5,5) gen zone brick size 5 5 5 p0 (0,0,5) edge 5.0 group Material ; Create Bin Zones
gen zone brick size 1 5 5 &
p0 (3,0,0) p1 add (3,0,0) p2 add (0,3,0) & p3 add (2,0,5) p4 add (3,6,0) p5 add (2,5,5) & p6 add (3,0,5) p7 add (3,6,5) gen zone brick size 1 5 5 &
p0 (5,0,5) p1 add (1,0,0) p2 add (0,5,0) & p3 add (0,0,5) p4 add (1,6,0) p5 add (0,5,5) & p6 add (1,0,5) p7 add (1,6,5) gen zone brick size 5 1 5 &
p0 (0,3,0) p1 add (3,0,0) p2 add (0,3,0) & p3 add (0,2,5) p4 add (6,3,0) p5 add (0,3,5) & p6 add (5,2,5) p7 add (6,3,5) gen zone brick size 5 1 5 &
p0 (0,5,5) p1 add (5,0,0) p2 add (0,1,0) & p3 add (0,0,5) p4 add (6,1,0) p5 add (0,1,5) & p6 add (5,0,5) p7 add (6,1,5) group Bin range group Material not ; Create named range synonyms range name=Bin group Bin
range name=Material group Material ; Assign models to groups model mohr range Material model elas range Bin
gen separate Material
interface 1 wrap Material Bin range plane ori 0 0 0 normal 1 -1 0 above
interface 2 wrap Material Bin range plane ori 0 0 0 normal 1 -1 0 below int 1 maxedge 0.55 int 2 maxedge 0.55
; Assign properties
prop shear 1e8 bulk 2e8 fric 30 range Material prop shear 1e8 bulk 2e8 range Bin ini den 2000
int 1 prop ks 2e9 kn 2e9 fric 15 int 2 prop ks 2e9 kn 2e9 fric 15 ; Assign Boundary Conditions
fix x range x -0.1 0.1 any x 5.9 6.1 any fix y range y -0.1 0.1 any y 5.9 6.1 any fix z range z -0.1 0.1 Bin
; Monitor histories hist unbal hist gp zdisp (6,6,10) hist gp zdisp (0,0,10) hist gp zdisp (0,0,0) ; Settings set large set grav 0,0,-10 ; Cycling step 4000 save bin.sav 9-5 ; ;------------------------------------------------------------------ rest pile_geom.sav model mohr range group clay prop bulk 8.333e7 shear 3.846e7 coh 30000 fric 0 range group clay model elas range group pile prop bulk 8.333e7 shear 3.846e7 range group pile interface 1 prop kn 1e8 ks 1e8 fric 20 coh 30000 interface 2 prop kn 1e8 ks 1e8 fric 20 coh 30000 ; ini dens 1230 range group clay ini dens 1230 range group pile model null range z -0.1 0.15 ; fix z range z -8.1 -7.9 fix x range x -8.1 -7.9 fix x range x 7.9 8.1 fix y range y -.1 .1 fix y range y 7.9 8.1 set grav 0 0 -10 ini szz 0. grad 0 0 12300. range z -5.5 0. ini szz 17600 grad 0 0 15500 range z -8 -5.5 ini sxx 0. grad 0 0 5271.4 range z -5.5 0. ini sxx 7542.86 grad 0 0 6642.86 range z -8 -5.5 ini sxx add 31428.6 grad 0 0 5714.3 range z -8 -5.5 ini syy 0. grad 0 0 5271.4 range z -5.5 0. ini syy 7542.86 grad 0 0 6642.86 range z -8 -5.5 ini syy add 31428.6 grad 0 0 5714.3 range z -8 -5.5 ; water density 1000 water table origin 0,0,-5.5 normal 0 0 -1 ini dens 1550 range z -8 -5.5 hist unbal ; solve rat 1.e-6 save pile0.sav ; model elas range group pile prop bulk 13.9e9 shear 10.4e9 range group pile ini dens 2500 range group pile call find_add.fis solve rat 1.e-6 save pile1.sav rest pile1.sav ;调用保存文件 ini state 0 ini xdis 0.0 ydis 0.0 zdis 0.0 ;位移清零
apply szz -0.4e6 range z 0.05 0.15 group pile ;桩顶加第一级荷载 solve save app0.4.sav print gp disp range id 1 ;输出第一级荷载下的桩顶位移，假定桩顶中心的id号为1 apply szz -0.6e6 range z 0.05 0.15 group pile ;桩顶加第二级荷载 solve save app0.6.sav print gp disp range id 1 ;输出第二级荷载下的桩顶位移 ???????????????? ;依次加载，直到桩破坏 ;-------------------------------------------------------------------- ;速度加载法 rest pile1.sav ini state 0 ini xdis 0 ydis 0 zdis 0 def zs_top ;检测桩顶竖向荷载 ad = top_head zftot = 0.0 loop while ad # null gp_pnt = mem(ad+1) zf = gp_zfunbal(gp_pnt)
zftot = zftot + zf ad = mem(ad) endloop
zs_top = zftot / 0.1414 end
fix z range z 0.05 .15 group pile ;固定桩顶速度，用速度来确定位移 def ramp
while_stepping if step < ncut then
udapp = float(step) * udmax / float(ncut) else
udapp = udmax endif
ad = top_head
loop while ad # null gp_pnt = mem(ad+1) gp_zvel(gp_pnt) = udapp ad = mem(ad) endloop end
hist gp zdis 0,0,0 hist gp zvel 0,0,0 hist zs_top
hist zone szz 0,0,-.1 set mech damp comb
set udmax = -1e-8 ncut 30000 step 225000 save pile2.sav ;
;--------------------------------------------------------------------
;位移控制法 def solve_steps loop n (1,21)
save_file = string(n) + '-step.sav' command step 40000 save save_file
pri zone stress ran id 2381 a id 2361 a id 2341 a ;输出桩顶网格单元的应力 endcommand endloop end
solve_steps
;----------------------------------------------------------------------------- 10-1
sel beam beg 0 0 0 end 2 0 0 nseg 2 sel beam beg 2 0 0 end 4 0 -1 nseg 3
sel beam id=2 beg 4 0 -1 end 5 0 -2 nseg 2 plot sel geo id on nod on scale 0.04 plot ad ax 10-2
sel node id=1 0 0 0 sel node id=2 2 0 0 sel node id=3 4 0 -1 sel node id=4 5 0 -2
sel beamsel id=1 cid=1 node 1 2 sel beamsel id=1 cid=2 node 2 3 sel beamsel id=1 cid=3 node 3 4 plot sel geo id on nod on scale 0.04 plot ad ax 10-3
def set_vals
ptA = 25.0 * sin( 40.0*degrad ) ptB = 25.0 * cos( 40.0*degrad ) end set_vals
gen zone cylinder p0=( 0.0, 0.0, 0.0 ) &
p1=( ptA, 0.0, ptB ) & p2=( 0.0, 25.0, 0.0 ) & p3=( 0.0, 0.0, 25.0 ) & p4=( ptA, 25.0, ptB ) & p5=( 0.0, 25.0, 25.0 ) & size=(1, 2, 2)
sel shell id=5 range cylinder end1=(0.0, 0.0,0.0) &
end2=(0.0,25.0,0.0) radius=24.5 not plot blo gro
plot ad sel geom black black cid on scale=0.03 plot ad ax pau
delete ; delete all zones
sel node init zpos add -25.0 10-4 new
gen zone brick size 6 8 8 model mohr
prop bulk 1e8 shear 0.3e8 fric 35 prop coh 1e10 tens 1e10 set grav 0 0 -9.81 ini dens 1000
fix x range x -0.1 0.1 fix x range x 5.9 6.1 fix y range y -0.1 0.1 fix y range y 7.9 8.1 fix z range z -0.1 0.1 hist n 5 hist unbal
set mech force 50 solve
save beam-brace0.sav ;
prop coh 1e3 tens 1e3
model null range x 2 4 y 2 6 z 5 10 set large
ini xdis 0 ydis 0 zdis 0
sel beam begin=( 2, 4, 8) end=( 4, 4, 8) nseg=2 sel beam prop emod=2.0e11 nu=0.30
sel beam prop XCArea=6e-3 XCIz=200e-6 XCIy=200e-6 XCJ=0.0 hist gp zdisp 4 4 8 solve
save beam-brace1.sav ;
plot create GravV
plot set plane dip 90 dd 0 origin 3 4 0 plot set rot 15 0 20
plot set center 2.5 4.2 4.0 plot set cap size 25
plot add cont disp plane behind shade on plot add sel beam force fx
plot add sel geom black black node=off shrinkfac=0.0 plot add axes plot show 10-5
;非全长锚固、预紧力锚杆（锚索）模拟
;方法1、通过删除－建立link链接来模拟托盘
gen zone radtun p0 0,0,0 p1 25,0,0 p2 0,50, 0 p3 0,0,25 size 4 25 4 10 dim 4 4 4 4 ratio 1 1 1 1.1 fill gen zone reflect normal 1 0 0 ori 0 0 0 gen zone reflect normal 0 0 1 ori 0 0 0 mo mohr
pro bulk 2.2e9 she 1.3e9 fric 30 coh 1.3e6 ten 1.5e5 ini dens 2000
fix x range x -25.1,-24.9 fix x range x 24.9 25.1 fix y range y 49.9 50.1 fix z range z -25.1 -24.9 fix z range z 24.9 25.1
sel cable id=1 beg 0, 0, 0 end 0 ,29, 0 nseg 10 sel cable id=1 beg 0,29,0 end 0,35,0 nseg 6
sel cable id=1 prop emod 2e10 ytension 310e3 xcarea 0.0004906 &
gr_coh 1 gr_k 1 gr_per 0.0785 range cid 1,10
sel cable id=1 prop emod 2e10 ytension 310e3 xcarea 0.0004906 &
gr_coh 10e5 gr_k 2e7 range cid 11,17 sel delete link range id 1 sel link id=100 1 target zone
sel link attach xdir=rigid ydir=rigid zdir=rigid xrdir=rigid yrdir=rigid zrdir=rigid range id 100 sel cable id=1 pretension 60e3 range cid 1,10
step 2000 sav 10-5.sav 10-6
;非全长锚固、预紧力锚杆（锚索）模拟 ;方法2、通过设置极大锚固剂参数模拟托盘
gen zone radtun p0 0,0,0 p1 25,0,0 p2 0,50, 0 p3 0,0,25 size 4 25 4 10 dim 4 4 4 4 ratio 1 1 1 1.1 fill gen zone reflect normal 1 0 0 ori 0 0 0 gen zone reflect normal 0 0 1 ori 0 0 0 mo mohr
pro bulk 2.2e9 she 1.3e9 fric 30 coh 1.3e6 ten 1.5e5 ini dens 2000
fix x range x -25.1,-24.9 fix x range x 24.9 25.1 fix y range y 49.9 50.1 fix z range z -25.1 -24.9 fix z range z 24.9 25.1
sel cable id=1 beg 0, 0, 0 end 0 ,29, 0 nseg 10 sel cable id=1 beg 0,29,0 end 0,35,0 nseg 6
sel cable prop emod 2e10 ytension 310e3 xcarea 0.0004906 &
gr_coh 1 gr_k 1 gr_per 0.0785 range cid 2,10
sel cable prop emod 2e10 ytension 310e3 xcarea 0.0004906 &
gr_coh 10e5 gr_k 2e7 range cid 11,17
sel cable prop emod 2e10 ytension 310e3 xcarea 0.0004906 &
gr_coh 10e8 gr_k 2e10 range cid 1,1
sel cable id=1 pretension 60e3 range cid 1,10 step 2000 sav 10-6.sav 10-7
;非全长锚固、预紧力锚杆（锚索）模拟
;方法3：借助别的结构单元（如liner单元）来模拟托盘
gen zone radtun p0 0,0,0 p1 25,0,0 p2 0,50, 0 p3 0,0,25 size 4 25 4 10 dim 4 4 4 4 ratio 1 1 1 1.1 fill gen zone reflect normal 1 0 0 ori 0 0 0 gen zone reflect normal 0 0 1 ori 0 0 0 mo mohr
pro bulk 2.2e9 she 1.3e9 fric 30 coh 1.3e6 ten 1.5e5 ini dens 2000
fix x range x -25.1,-24.9 fix x range x 24.9 25.1 fix y range y 49.9 50.1 fix z range z -25.1 -24.9 fix z range z 24.9 25.1
sel cable id=1 beg 0, 0, 0 end 0 ,29, 0 nseg 10 sel cable id=1 beg 0,29,0 end 0,35,0 nseg 6 sel cable id=1 prop emod 2e10 ytension 310e3 xcarea 0.0004906 &
gr_coh 1 gr_k 1 gr_per 0.0785 range cid 1,10
sel cable id=1 prop emod 2e10 ytension 310e3 xcarea 0.0004906 &
gr_coh 10e5 gr_k 2e7 range cid 11,17 sel liner range y=-.1, .1 x=-1,1 z=-1,1
sel liner PROP iso=( 25e9, 0.15) thick=0.1 ; concrete sel liner PROP cs_nk=8e8 cs_sk=8e8 &
cs_ncut=0.0 cs_scoh=0.0 cs_scohres=0.0 cs_sfric=0.0 sel delete link range id 1
sel link id=100 1 target node tgt_num 18
sel link attach xdir=rigid ydir=rigid zdir=rigid xrdir=rigid yrdir=rigid zrdir=rigid range id 100 sel cable id=1 pretension 60e3 range cid 1,10
step 2000 sav 10-7.sav 10-8 n
title Structure_dynamic_analysis_lakewater config dyn
sel pile id=1 beg 0 0 0 end 0 0 1 sel pile prop dens 2400 &
Emod 1.0e10 Nu 0.3 XCArea 0.3 &
XCJ 0.16375 XCIy 0.00625 XCIz 0.01575 & Per 2.8 &
CS_sK 1.3e11 CS_sCoh 0.0 CS_sFric 10.0 &
CS_nK 1.3e11 CS_nCoh 0.0 CS_nFric 0.0 CS_nGap off
def f1
whilestepping
f0=10000*sin(10*dytime) np = nd_head
loop while np # null if nd_pos(np,1,3)=1 nd_apply(np,1)=f0 endif
np = nd_next(np) endloop end
sel node fix x y z xr yr zr ran id=1 sel set damp combined
plo cre pile plo current pile
plo set back black fore white mag 0.8 plo add sel geo id on sca .04 magf 1e3 plo add sel fapp lgreen magf 1e3
plo add sel pile mom my lblue lred magf 1e3 axe yel set movie avi step 100 file pile.avi movie start sol age 1 movie finish 11-1 new
conf dyn ;打开动力计算功能 gen zone brick size 10 5 10 mod elas
mod null range x=0,5 z=5,10 ;删除部分网格 fix z range x=-.1 .1 z=.1 10.1 ;设置静力边界条件 fix z range x=9.9,10.1 z=.1 10.1 fix y range y=-.1 .1 fix y range y=4.9 5.1
prop bulk 2e8 shear 1e8 ;设置土体参数 prop bulk 4e9 shear 2e9 range x=5,6 z=5,10 ;设置墙体参数（土体参数的20倍）
ini dens 2000 ;设置密度
def setup ;动荷载中的变量赋值 freq = 1.0
omega = 2.0 * pi * freq old_time = clock end
setup ;执行变量赋值 def wave ;定义动荷载函数 wave = sin(omega * dytime) ;定义动荷载变量 end
apply xvel = 1 hist wave range z=-.1 .1 ;施加动荷载 apply zvel = 0 range z=-.1 .1 hist gp xvel 5,2,0 hist gp xvel 5,2,10 hist gp zvel 5,2,10 hist dytime
def tim ;估算程序运行的时间 tim = 0.01 * (clock - old_time) end set dyn multi on ;设置动态多步 solve age 1.0
print tim ;输出计算时间
print dyn ;输出动力计算相关信息 save mult1.sav 11-2 new
config dyn
gen zone brick size 1,1,50 model elas
prop shear 1e7 bulk 2e7 ini dens 1000 def setup
omega = 2.0 * pi * freq pulse = 1.0 / freq end
set freq=4.0 setup def wave
if dytime > pulse wave = 0.0 else
wave = 0.5 * (1.0 - cos(omega * dytime)) endif end
range name bottom z=-.1 .1
fix z range z=.5 55 ;将上部网格都施加数值向约束 apply dquiet squiet range bottom
apply sxz -2e5 hist wave syz 0.0 szz 0.0 range bottom ;-2e5的系数来源于 的值
apply nvel 0 plane norm 0,0,1 range bottom hist gp xvel 0,0,0 hist gp xvel 0,0,25 hist gp xvel 0,0,50 hist dytime hist wave
plot create hhh
plot add hist 1 2 3 vs 4 plot show solve age 2
11-3 new
;第一步：静力计算阶段 config dyn set dyn off
gen zone brick size 6 3 2
gen zone brick size 2 3 2 p0 0 0 2 gen zone brick size 2 3 2 p0 4 0 2 gen zone wedge size 1 3 2 p0 2 0 2
gen zone wedge size 1 3 2 p0 4 3 2 p1 3 3 2 p2 4 0 2 p3 4 3 4 &
p4 3 0 2 p5 4 0 4
model elastic
prop bulk 66667 shear 40000 ini dens 0.0025 set grav 0 0 -10
fix x range x -0.01 0.01 fix x range x 5.99 6.01 fix y range y -0.01 0.01 fix y range y 2.99 3.01 fix z range z -0.1 0.1 hist unbal solve
save 11-3_1.sav
;第二步：动力计算阶段 set dyn on def iniwave
per = 0.01 end iniwave def wave
wave = 0.5 * (1.0 - cos (2*pi*dytime/per)) end
free x y z ran z -0.1 0.1 ;去掉模型底部原有的静力条件
apply nquiet squiet dquiet ran z -0.1 0.1 ;静态边界条件
apply dstress 1.0 hist wave ran z -0.1 0.1 ;加动力荷载 apply ff ;施加自由场边界条件 group ff_corner
group ff_side ran x 0 6 group ff_side ran y 0 3
group main_grid ran x 0 6 y 0 3
set dyn time = 0 ;设置动力计算从0s开始 hist reset ;清空已有的历史信息 hist unbal
hist dytime ; 主体网格 hist gp xvel 2 1 0
hist gp xvel 2 1 5.0 ; 柱体网格 hist gp xvel -1 -1 0
hist gp xvel -1 -1 5 ;平行于y方向的二维自由场网格 hist gp xvel -1 0 0 hist gp xvel -1 0 5.0
;平行于x方向的二维自由场网格 hist gp xvel 2 -1 0 hist gp xvel 2 -1 5.0 solve age 0.015 save 11-3_2.sav 11-4 conf dy
gen zone brick size 3,3,3 model elas
prop bulk 1e8 shear 0.3e8 ini dens 1000 fix z range z -.1,.1
set dyn=on, grav 0 0 -10, hist_rep=1 hist gp zdisp 3.0,1.5,3.0 hist dytime plot create hh plot add his 1 vs 2 save 11-4_1.sav
;保存文件，在后续计算中需调用该文件 cyc 150 11-5
;（1）质量分量和刚度分量共同作用 rest 11-4_1.sav
set dyn damp rayleigh 1 22.8 solve age=0.2 title
vertical displacement versus time (mass & stiffness damping) plot show pause
;（2）只有质量分量 rest 11-4_1.sav
set dyn damp rayleigh 2 22.8 mass solve age=0.08 title
vertical displacement versus time (mass damping only) plot show pause
;（3）只有刚度分量 rest 11-4_1.sav
set dyn damp rayleigh 2 22.8 stiffness solve age=0.08 title
vertical displacement versus time (stiffness damping only) plot show 11-6
rest 11-4_1.sav
set dyn damp rayleigh 0.05 22.8 set hist_rep=5 solve age=0.5 title
vertical displacement versus time (5% Rayleigh damping) plot show pause
rest 11-4_1.sav
set dyn damp local 0.1571 ; = pi * 0.05 set hist_rep=5 solve age=0.5 title
vertical displacement versus time (5% Local damping) plot show 11-7
;振动台试验的例子 new
config dynamic fluid def model_dim h_R = 0
h_R1 = h_R + 1.0 end
model_dim gen zon bri p0 0 0 -10 p1 30 0 -10 p2 0 1 -10 p3 0 0 0 p4 30 1 -10 p5 0 1 0 p6 30 0 h_R p7 30 1 h_R size 30 1 10 group sand
gen zon bri p0 0 0 0 p1 30 0 h_R p2 0 1 0 p3 0 0 1 p4 30 1 h_R p5 0 1 1 p6 30 0 h_R1 p7 30 1 h_R1 size 30 1 1 group top
;gen zon bri p0 0 0 -.5 p1 3 0 -.5 p2 0 1 -.5 p3 0 0 0 p4 3 1 -.5 p5 0 1 0 p6 3 0 0 p7 3 1 0 size 30 1 10 model elastic
prop bulk=3e7 shear=1e7 fric=35 ini dens 2000 model fl_iso
prop poro 0.5 perm 1e-8 ini fmod 2e8 fdens 1000
ini pp 0 grad 0 0 -10e3 ran z 0 -10.0 fix z ran z -9.9 -10.1 fix x ran x -.1 .1 fix x ran x 29.9 30.1 fix y
set grav 10
set fluid off dyn off ini fmod 0
set mech rat 1e-6 solve
def ini_conf _k0 = 1.0
pnt = zone_head loop while pnt # null
val = _k0 * z_szz(pnt) + (_k0-1.) * z_pp(pnt) z_sxx(pnt) = val z_syy(pnt) = val pnt=z_next(pnt) end_loop end ini_conf solve
save 11-7.sav 11-8
rest 11-7.sav
set dyn on fluid on ini fmod 2e8 set fluid pcut on
model finn ran gro sand
prop bulk=3e7 shear=1e7 fric=35 ran gro sand endloop ini dens 2000 ran gro sand end
prop ff_latency=50 ff_switch = 0 ff_c1=0.8 solve_ages ff_c2=0.79 ff_c3=0.45 ff_c4=0.73 ran gro sand ;扭剪试验结果
def setup freq=5.0 ampl=2 omega = 2.0 * pi * freq end setup def sine_wave vv = 9.36e-2 * sin(omega*dytime) if dytime < 2.0 sine_wave = dytime / 2.0 * vv else if dytime < 20.0 sine_wave = vv else if dytime <= 30. sine_wave = (30.0 - dytime) / 10.0 * vv endif endif endif if dytime > 30.0 sine_wave = 0.0 endif end free x apply xvel = 1.0 hist sine_wave ran z -9.9 -10.1 apply xvel = 1.0 hist sine_wave ran x -.1 .1 apply xvel = 1.0 hist sine_wave ran x 29.9 30.1 set dyn damp local .314 call ppr.dat set hist_rep 100 set large ;set dyn dt 3e-4 set mech rat 1e-20 def solve_ages loop n(1,39) save_file = '11-8_'+string(n)+'s.sav' command sol age n ;save save_file endcommand save save_file
hist write 20 21 22 23 24 30 31 32 33 34 vs 2 file 10-8_Outfile_pp.txt hist write 40 41 42 43 44 vs 2 file 10-8_Outfile_xdis.txt 12-1 gen zon bri size 3 3 3 ini pp 30e3 grad 0 0 -10e3 plot con pp outline on ;outline的作用是在云图中显示网格的轮廓 12-2 gen zon bri size 3 3 3 water density 1000 set grav 10 water table ori 0 0 3 norm 0 0 1 12-3 gen zon bri size 3 3 3 ini pp 30e3 grad 0 0 -10e3 ran x -.1 .1 fix pp ran x -.1 .1 plot con pp outline on： 显示节点上的孔压云图 plot bcon pp outline on： 显示单元中心点处平均孔压块图 12-4 gen zon bri size 3 3 3 app pp 30e3 grad 0 0 -10e3 ran x -.1 .1 plot con pp outline on： 显示节点上的孔压云图 plot bcon pp outline on： 显示单元中心点处平均孔压块图 12-5 config fluid gen zone brick size 20 1 10 model mohr prop bulk 8.33e6 shear 3.85e6 fric 15 coh 10e3 tens 1e10; fix x range x -.1 .1 fix x range x 19.9 20.1 fix x y z range z -.1 .1
fix y
; --- apply load slowly --- def ramp
ramp = min(1.0,float(step)/200.0) end
apply nstress = -40e3 hist ramp range x -.1 3.1 z 9.9 10.1 ; --- fluid flow model --- model fl_iso ini fmod 2e9
; --- pore pressure fixed at zero at the surface --- fix pp 0 range z 9.9 10.1 ; --- settings --- set fl off
; --- histories --- hist gp pp 2,.5,9 ; --- test --- step 750 save load.sav
plot set plane ori 0 0.5 0 norm 0 1 0 plot con pp plane ou on plot add fap red plane 12-6
config fluid set fluid off
gen zon brick p0 0 0 -10 size 20 1 10
gen zon brick p0 5 0 0 p1 15 0 0 p2 5 1 0 p3 9 0 5 p4 15 1 0 p5 9 1 5 p6 11 0 5 p7 11 1 5 & size 10 1 5 group soil
group dam ran x 5 7 z -5 0
group dam ran id 201 a id 211 a id 221 a id 231 a id 241 a
group dam ran id 202 a id 212 a id 222 a id 232 a id 242 a m e
prop bu 3e7 sh 1e7
ini pp 0 grad 0 0 -10e3 ran z 0 -10 ini dens 2000 model fl_iso
prop por 0.5 perm 1e-10 ini fden 1000 ften -1e10 ini sat 0.0 ran z 0 5
model fl_null ran gro dam ;ini pp 0 ran gro dam fix z ran z -10 fix x ran x 0 fix x ran x 20 fix y
set grav 10 solve
save elastic.sav 12-7
rest elastic.sav
ini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0
app nstress -40e3 grad 0 0 10e3 ran z 0 4 x 0 9 solve
save pressure.sav 12-8
rest pressure.sav set fluid on mech off
ini fmod 2e3 ften 0.0 ran gro soil ini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0
app pp 40e3 grad 0 0 -10e3 ran z 0 4 x 0 9 app pp 0 ran z 0 x 15 20 hist id=10 zone pp id 215 solve 12-9
config fluid set fluid biot on
gen zone bri p0 0 0 0 p1 5 0 0 p2 0 2 0 p3 0 0 20 size 3 1 10
group Face range z 18 20
group Soft range z 10 18 ;soft soil 8m group Clay range z 0 10 model mohr
prop bulk 33e6 shear 7e6 c 0 f 30 range group Face
prop bulk .83e6 shear .17e6 c 7.0e3 f 5 range group Soft
prop bulk 6.66e6 shear 1.42e6 c 8.0e3 f 20 range group Clay ini dens 1500 model fl_iso
ini ftens -5e5 ini fdens 1000 ini biot_mod 4e9
prop perm 10e-13 biot_c 1 prop perm 10e-10 ran grou Face fix x y z range z -0.1 0.1 fix y rang y -.01 0.01 fix y rang y 1.9 2.1 fix x range x -.01 .01 fix x range x 4.9 5.1 ;Initial Stress K0=1 set grav 0 0 -10
ini szz -400e3 grad 0 0 20e3 ini syy -400e3 grad 0 0 20e3 ini sxx -400e3 grad 0 0 20e3 ini pp 200e3 grad 0 0 -10e3 set fluid off step 5000
ini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0 ;for fixed pp on the PVD position def fixPP loop i(1,10)
pv = -80e3 +(i-1)*10e3 zt = 20.1-(i-1) zb = 20.1-i command
fix pp pv rang x -.1 .1 z zb zt end_command end_loop end fixpp step 5000 set fluid on
ini biot_mod 4.3e7 rang group Soft not ;modify biot_mod
ini biot_mod 9.0e6 rang group Soft set mech force 1.5e3
set mech subs 100 auto ;slave set fluid subs 30 ; --- histories --- def day
day=fltime/24/3600 end day
hist id=1 day hist id=2 gp zd 2.5,0,20 hist id=3 gp zd 2.5,0,10 sol age 7.776e6 13-1
config cppudm ;必须设置cppudm的选项才可以进行模型载入
model load DuncanChang.dll
gen zon bri p1 .6 0 0 p2 0 .6 0 p3 0 0 .6 size 1 1 1;三轴试验尺寸：0.6*0.6*0.6 m model duncan
prop cohesion 110e3 friction 48.5 fricDel 0.0 ratiofail 0.79 ke 704 ne 0.38 kb 303 mb 0.18 kur 844.8
fix z ran z -.01 .01 ;模型底部边界的竖直向速度约束 def sigma3 ;定义一个sigma变量便于进行多工况计算 sigma3=-600e3 end sigma3
;设置初始应力
app nstress sigma3 ran x -.01 .01 app nstress sigma3 ran x .59 .61 app nstress sigma3 ran y -.01 .01 app nstress sigma3 ran y .59 .61 app nstress sigma3 ran z .59 .61
;ini den 2190 ;为了与理论值作比较，不考虑重力、密度的作用 ini szz sigma3 ini syy sigma3 ini sxx sigma3
solve ;得到加载前的初始应力
ini xdis 0 ydis 0 zdis 0 xvel 0 yvel 0 zvel 0 ;位移和速度清零
tab 1 name loads ;建立一个空表，用来保存荷载-沉降曲线 plo tab 1
;第1次加卸载 def load1
p_gp = gp_near(0,0,0.6) loop n(1,50)
zss_load= sigma3 - float(n)*20e3 ;加载1000kPa command
app nstress zss_load ran z .59 .61 solve
endcommand
z_dis = -1.0 * gp_zdisp(p_gp) / 0.6 ;得到应变 z_load = (sigma3 - zss_load) ;得到主应力差
command
tab 1 z_dis z_load ;分别保存应变和主应力差 endcommand end_loop end load1
save 13-1.sav 13-2
rest 13-1.sav def unload1
p_gp = gp_near(0,0,0.6) loop m(1,25)
zss_load= sigma3 - 1000e3 + float(m)*2e3 ;卸载500kPa
command
app nstress zss_load ran z .59 .61 solve
endcommand
z_dis = -1.0 * gp_zdisp(p_gp) / 0.6 ;得到应变 z_load = (sigma3 - zss_load) ;得到主应力差
command
tab 1 z_dis z_load ;分别保存应变和主应力差 endcommand endloop end unload1 pau
set log on
set logfile 1.log pri tab 1 set log off
;第2次加卸载 def load2
p_gp = gp_near(0,0,0.6) loop n(1,75)
zss_load= sigma3 - 1000e3 + 500e3 - float(n)*20e3 ;加载至2000kPa command
app nstress zss_load ran z .59 .61 solve
endcommand
z_dis = -1.0 * gp_zdisp(p_gp) / 0.6 ;得到应变 z_load = (sigma3 - zss_load) ;得到主应力差
command
tab 1 z_dis z_load ;分别保存应变和主应力差 endcommand end_loop end load2
save 13-2.sav 14-1 n
;=================================== ;建立网格模型 gen zone brick &
p0 0 0 0 p1 2 0 0 p2 0 0.5 0 p3 0 0 3 size 3 1 3 gen zone brick &
p0 2 0 0 p1 20 0 0 p2 2 0.5 0 p3 2 0 3 & size 17 1 3 ratio 1.03 1 1 gen zone brick &
p0 2 0 3 p1 20 0 3 p2 2 0.5 3 p3 12 0 13 & p4 20 0.5 3 p5 12 0.5 13 p6 20 0 13 & p7 20 0.5 13 size 17 1 17 ratio 1.03 1 1
;==================================== ;设置边界条件
fix x y z range z -0.1 0.1 fix x range x 19.9 20.1 fix x range x -0.1 0.1 fix y
;==================================== ;初始地应力的生成 model elas
prop density 2000 bulk 3e9 shear 1e9 set gravity 0 0 -10.0 solve
ini xdisp 0 ydisp 0 zdisp 0 ini xvel 0 yvel 0 zvel 0
;===================================
;安全系数求解 model mohr
prop density 2000 bulk 1.0e8 shear 3.0e7 & coh 12380.0 tens 1.0e6 fric 20 dila 20 solve fos file slope3dfos1.sav associated 14-2 n
;===================================== ;建立网格模型 gen zone brick &
p0 0 0 0 p1 2 0 0 p2 0 0.5 0 p3 0 0 3 size 3 1 3 gen zone brick &
p0 2 0 0 p1 20 0 0 p2 2 0.5 0 p3 2 0 3 size 17 1 3 & ratio 1.03 1 1
gen zone brick &
p0 2 0 3 p1 20 0 3 p2 2 0.5 3 p3 12 0 13 & p4 20 0.5 3 p5 12 0.5 13 p6 20 0 13 & p7 20 0.5 13 size 17 1 17 ratio 1.03 1 1
;**********************************************
;自定义强度折减法 def SSR
;===================================== ;定义有关参数及循环终止条件 ait1=0.02 k11=1.0 k12=2.0 ks=(k11+k12)/2
loop while (k12-k11)>ait1 coh1=12380/ks
fri1=(atan((tan(20*pi/180))/ks))*180/pi dila1=20.0 ten1=1e6 grav0=-10
dens1=2000 K1=1e8 G1=3e7
;===================================== ;折减的实现过程 command model null
;初始应力场的生成 model elastic pro bulk 1e10 she 3e9 dens dens1 fix x y z range z -0.1 0.1 fix x range x 19.9 20.1 fix x range x -0.1 0.1 fix y
set grav 0 0 grav0 solve
ini xdisp 0 ydisp 0 zdisp 0 ini xvel 0 yvel 0 zvel 0 ;塑性阶段求解 model mohr
pro bulk K1 she G1 dens dens1 coh coh1 & friction fri1 dil dila1 tens ten1 set mech ratio 9.8e-6 solve step 10000 endcommand
;二分法的实现过程 if mech_ratio<1.0e-5 k11=ks k12=k12 else
k12=ks k11=k11 endif
ks=(k11+k12)/2 endloop
;===================================== ;计算结果的保存 fosfile0='_fos'+'.sav' command
save fosfile0 endcommand end
;**********************************************
;程序执行及结果显示 SSR
pr ks 15-1 n
;导入分组网格 impgrid 1
group soil range group rock not
;定义材料性质 model mohr
pro den 2200 bulk 4e7 she 1.5e7 co 8e4 & fric 28 ten 9e3 dil 10 range group soil
pro den 2900 bulk 2.5e10 she 1.2e10 co 4e6 & fric 45 ten 2e6 dil 15 range group rock ;设置边界条件 fix x y z range z -0.1 0.1 fix x y z range z 7.9 8.1 fix y ;加载
app nstress -5e5 range x -0.1 0.1 app nstress -5e5 range x 4.9 5.1 ;设置初始条件
ini zvel -2.5e-5 range z 7.9 8.1 ;定义轴向应力求解函数 def axi_stress f_accum=0.0 pnt = gp_head
loop while pnt # null
if gp_zpos(pnt) <0.15 then
f_accum = f_accum + gp_zfunbal(pnt) endif
pnt = gp_next(pnt) endloop
axi_stress = f_accum / 0.5 end
;定义侧向应力求解函数 def lat_stress g_accum=0.0 pnt = gp_head
loop while pnt # null
if gp_xpos(pnt) <0.1 then
g_accum = g_accum + gp_xfunbal(pnt) endif
pnt = gp_next(pnt) endloop
lat_stress = g_accum / 0.8 end
;跟踪记录有关重要变量 hist n 1
hist axi_stress hist lat_stress hist gp zdis 0 0 8 hist gp xdis 0 0 4 hist gp xdis 0 5 4 hist unbal
;绘制应力-位移曲线
plot hist -1 vs -3 ;axial stress vs axial displace ;求解
step 15000
;结果保存及数据输出 save 1.sav
hist write 1 file axi_stress.dat hist write 2 file lat_stress.dat hist write 3 file axi_dis.dat hist write 4 file lat_dis1.dat hist write 5 file lat_dis2.dat hist write 6 file unbalance.dat 15-2 n
;=================================== ;导入网格数据并生成“辅助”界面单元 impgrid ww gro water
interface 1 face ran x -4.9 604.9 y 0.1 300 z -4.9 704.9 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 aa
group soil range group 2 any group 4 any group rock range group 1 any group 3 any ;生成潜水面 water_table
;删除为“辅助”网格单元和界面单元 int 1 dele
dele range group water
;初始化材料密度: 饱和的和非饱和的 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
;施加边界约束条件 fix x y z rang y -0.1 0.1
fix z range x 0 600 z -0.1 0.1 fix x range x -0.1 0.1 z 0 700
fix z range x 0 600 z 699.9 700.1 fix x range x 599.9 600.1 z 0 700
;通过弹性求解生成初始应力场 model elas
pro bulk 1e10 she 1e10 range group soil pro bulk 1e10 she 1e10 range group rock solve fo 10
;位移和速度归零
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
set mech ratio 1.0e-6 solve