(完整word版)用flac3d模拟基坑开挖

new
;网格建立
;;;;;;;;;;;;;;;;;;;;;;;;;
gen zone brick p0 90 0 -30 p1 202 0 -30 p2 90 4 -30 p3 90 0 0 size 112 4 30 ratio 1 1 1
gen zone brick p0 90 0 -30 p1 90 0 0 p2 90 4 -30 p3 0 0 -30 size 30 4 25 ratio 1 1 1.1
gen zone brick p0 90 0 -30 p1 0 0 -30 p2 90 4 -30 p3 90 0 -75 size 25 4 18 ratio 1.1 1 1.1
gen zone brick p0 90 0 -30 p1 90 0 -75 p2 90 4 -30 p3 202 0 -30 size 18 4 112 ratio 1.1 1 1 gen zone brick p0 202 0 -30 p1 292 0 -30 p2 202 4 -30 p3 202 0 0 size 25 4 30 ratio 1.1 1 1 gen zone brick p0 202 0 -30 p1 202 0 -75 p2 202 4 -30 p3 292 0 -30 size 18 4 25 ratio 1.1 1 1.1
;分组
;;;;;;;;;;;;;;;;;;;;;;;;;;
group 1 range x 90 110 y 0 4 z -30 0
group 1 range x 180 202 y 0 4 z -30 0
group 2 range group 1 not
;建立连续墙单元
;;;;;;;;;;;;;;;;;;;;;;;;;;
gen separate 1
gen merge 1e-4 range x 90 110 y 0 4 z -30.1 -29.9
gen merge 1e-4 range x 180 202 y 0 4 z -30.1 -29.9
attach face range x 89.99 90.01 y 0.0 4.0 z -29.9 0
attach face range x 109.99 110.01 y 0.0 4.0 z -29.9 0
attach face range x 179.99 180.01 y 0.0 4.0 z -29.9 0
attach face range x 201.99 202.01 y 0.0 4.0 z -29.9 0
sel liner id 1 crossdiag group 2 range x 89.9 90.1 y -0.1 4.1 z -30.1 0.1
sel liner id 2 crossdiag group 2 range x 109.9 110.1 y -0.1 4.1 z -30.1 0.1
sel liner id 3 crossdiag group 2 range x 179.9 180.1 y -0.1 4.1 z -30.1 0.1
sel liner id 4 crossdiag group 2 range x 201.9 202.1 y -0.1 4.1 z -30.1 0.1
sel liner id 1 prop isotropic (2.0e10, 0.20) thickness 0.8 density 2.5e3 &
cs_nk=4e9 cs_sk=4e9 &
cs_ncut=4e7 cs_scoh=4e7 cs_scohres=0 cs_sfric=20.0 &
range x 89.9 90.1 y -0.1 4.1 z -30.1 0.1
sel liner id 2 prop isotropic (2.0e10, 0.20) thickness 0.8 density 2.5e3 &
cs_nk=4e9 cs_sk=4e9 &
cs_ncut=4e7 cs_scoh=4e7 cs_scohres=0 cs_sfric=20.0 &
range x 109.9 110.1 y -0.1 4.1 z -30.1 0.1
sel liner id 3 prop isotropic (2.0e10, 0.20) thickness 0.8 density 2.5e3 &
cs_nk=4e9 cs_sk=4e9 &
cs_ncut=4e7 cs_scoh=4e7 cs_scohres=0 cs_sfric=20.0 &
range x 179.9 180.1 y -0.1 4.1 z -30.1 0.1
sel liner id 4 prop isotropic (2.0e10, 0.20) thickness 0.8 density 2.5e3 &
cs_nk=4e9 cs_sk=4e9 &
cs_ncut=4e7 cs_scoh=4e7 cs_scohres=0 cs_sfric=20.0 &
range x 201.9 202.1 y -0.1 4.1 z -30.1 0.1
;定义支撑结构
;;;;;;;;;;;;;;;;;;;;;;;
def struct_install1
loop i(1,3)
structx_zz=-1.0*5.0*(i-1)
structx_xx0=90.0
structx_xx1=110.0
structx_yy=2.0
command
sel beam id=2 begin (structx_xx0,structx_yy,structx_zz) end (structx_xx1,structx_yy,structx_zz) nseg=10
sel beam id=2 prop dens=0.000 emod=1.0e-6 nu=0.2 &
xcarea=0.80 xcj=10.94e-2 xciy=6.67e-2 xciz=4.27e-2 ydirection=(0 0 -1) ;1000x800
endcommand
endloop
end
struct_install1
def struct_install2
loop i(1,3)
structx_zz=-1.0*5.0*(i-1)
structx_xx0=180.0
structx_xx1=202.0
structx_yy=2.0
command
sel beam id=3 begin (structx_xx0,structx_yy,structx_zz) end (structx_xx1,structx_yy,structx_zz) nseg=11
sel beam id=3 prop dens=0.000 emod=1.0e-6 nu=0.2 &
xcarea=0.80 xcj=10.94e-2 xciy=6.67e-2 xciz=4.27e-2 ydirection=(0 0 -1) ;1000x800
endcommand
endloop
end
struct_install2
;建立结构单元分组
;;;;;;;;;;;;;;;;;;;;;;;;;;;
sel group linerwall range sel liner
sel group struct1 range sel beam x (90.0 110.0) z (-0.1 0.1)
sel group struct2 range sel beam x (90.0 110.0) z (-5.1 -4.9)
sel group struct3 range sel beam x (90.0 110.0) z (-10.1 -9.9)
sel group struct4 range sel beam x (180.0 202.0) z (-0.1 0.1)
sel group struct5 range sel beam x (180.0 202.0) z (-5.1 -4.9)
sel group struct6 range sel beam x (180.0 202.0) z (-10.1 -9.9)
;删除beam单元的link
sel dele link range sel beam z (-30 0)
;建立liner间的节点间的刚性link
def merge_link0
node_num=0
node_pnt0 = nd_head
loop while node_pnt0 # null ;寻找总节点数，注：不能自己任生成node，程序缺省的方式为连续生成无不连续
node_num = node_num+1
node_pnt0 = nd_next(node_pnt0)
endloop
node_num_minus1 = node_num-1
link_id=30000
loop ii (1,node_num_minus1)
node_pnt1 = nd_find(ii)
xxa = nd_pos(node_pnt1,2,1)
yya = nd_pos(node_pnt1,2,2)
zza = nd_pos(node_pnt1,2,3)
ii_plus1 = ii+1
loop jj (ii_plus1,node_num)
node_pnt2 = nd_find(jj)
xxb = nd_pos(node_pnt2,2,1)
yyb = nd_pos(node_pnt2,2,2)
zzb = nd_pos(node_pnt2,2,3)
node_dist = sqrt((xxa-xxb)^2+(yya-yyb)^2+(zza-zzb)^2)
dist_tol = 1e-1
if node_dist <= dist_tol then
link_pnt1 = nd_link(node_pnt1)
link_pnt2 = nd_link(node_pnt2)
;if link_pnt1 # null then
; temp1 = lk_delete(link_pnt1)
;endif
if link_pnt2 # null then
temp2 = lk_delete(link_pnt2)
endif
link_id = link_id+1
command ;生成新link(6自由度全固结),大的node的id作为target node,小的node的id作为source node，需注意不同情况下的灵活调整
sel set link node_tol=dist_tol
sel link id=link_id jj target = node tgt_num =ii ;指定link的ID
;sel link ii target = node tgt_num = jj ;不指定link的id，自动生成
sel link attach xdir=rigid ydir=rigid zdir=rigid xrdir=rigid yrdir=rigid zrdir=rigid range id=link_id
endcommand
endif
endloop
endloop
end
merge_link0
;设置土层材料参数
;;;;;;;;;;;;;;;;;;;;;;;;;;;;
def b_s_mod
b_mod =e_mod/(3.0*(1.0-2.0*p_ratio))
s_mod =e_mod/(2.0*(1.0+p_ratio))
end
model elastic
set e_mod 100e6
set p_ratio 0.3
b_s_mod
prop bu=b_mod sh=s_mod
ini dens 1800 range z -75 0
def ini_szz
szz0=0
szzgrad=1800*10
command
ini szz add szz0 grad 0 0 szzgrad range z -75 0
endcommand
end
ini_szz
def ini_sxx_syy
pnt=zone_head
loop while pnt # null
val=k0*z_szz(pnt)
z_sxx(pnt)=val
z_syy(pnt)=val
pnt=z_next(pnt)
endloop
end
set k0=0.50
ini_sxx_syy
;定义边界处的结构边界条件;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
cyc 0
sel node local xdir=(0,1,0) ydir=(0,0,1) range x 89.9 90.1 y -0.1 4.1 z -30.1 0.1 sel node local xdir=(0,1,0) ydir=(0,0,-1) range x 109.9 110.1 y -0.1 4.1 z -30.1 0.1 sel node local xdir=(0,1,0) ydir=(0,0,1) range x 179.9 180.1 y -0.1 4.1 z -30.1 0.1 sel node local xdir=(0,1,0) ydir=(0,0,-1) range x 201.9 202.1 y -0.1 4.1 z -30.1 0.1
sel node fix lsys range x 89.9 90.1 y -0.1 0.1 z -30.1 0.1
sel node fix lsys range x 89.9 90.1 y 3.9 4.1 z -30.1 0.1
sel node fix lsys range x 109.9 110.1 y -0.1 0.1 z -30.1 0.1
sel node fix lsys range x 109.9 110.1 y 3.9 4.1 z -30.1 0.1
sel node fix lsys range x 179.9 180.1 y -0.1 0.1 z -30.1 0.1
sel node fix lsys range x 179.9 180.1 y 3.9 4.1 z -30.1 0.1
sel node fix lsys range x 201.9 202.1 y -0.1 0.1 z -30.1 0.1
sel node fix lsys range x 201.9 202.1 y 3.9 4.1 z -30.1 0.1
sel node fix x yr zr range x 89.9 90.1 y -0.1 0.1 z -30.1 0.1
sel node fix x yr zr range x 89.9 90.1 y 3.9 4.1 z -30.1 0.1
sel node fix x yr zr range x 109.9 110.1 y -0.1 0.1 z -30.1 0.1
sel node fix x yr zr range x 109.9 110.1 y 3.9 4.1 z -30.1 0.1
sel node fix x yr zr range x 179.9 180.1 y -0.1 0.1 z -30.1 0.1
sel node fix x yr zr range x 179.9 180.1 y 3.9 4.1 z -30.1 0.1
sel node fix x yr zr range x 201.9 202.1 y -0.1 0.1 z -30.1 0.1
sel node fix x yr zr range x 201.9 202.1 y 3.9 4.1 z -30.1 0.1
sel node fix y range x 89.9 90.1 y 0.0 4.0 z -0.1 0.1
sel node fix y range x 109.9 110.1 y 0.0 4.0 z -0.1 0.1
sel node fix y range x 179.9 180.1 y 0.0 4.0 z -0.1 0.1
sel node fix y range x 201.9 202.1 y 0.0 4.0 z -0.1 0.1
;set plot meta
;plot set rot 20 0 30 ba wh color=on cent=(10 20 0) mag=3.81
;set outp node_local_sys.wmf
;plot add sel geom black red link=off node=off id=off shrink=0 scale=0.03 nodesys=on range group linerwall any group struct1 any
;pl ha
;固定边界条件
;;;;;;;;;;;;;;;;;;;;;;;;;;
fix x range x -0.1 0.1
fix x range x 291.9 292.1
fix y range y -0.1 0.1
fix y range y 3.9 4.1
fix x y z range z -75.1 -74.9
set grav 0,0,-10
solve
save elas.sav
;删除侧面内外土体间的连接约束;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
attach delete range x 89.99 90.01 y 0.0 4.0 z -29.9 0
attach delete range x 109.99 110.01 y 0.0 4.0 z -29.9 0
attach delete range x 179.99 180.01 y 0.0 4.0 z -29.9 0
attach delete range x 201.99 202.01 y 0.0 4.0 z -29.9 0
;在墙内土体的外侧建立接触面;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
interface 1 face range group 1 x 89.99 90.01 y 0.0 4.0 z -29.9 0
interface 2 face range group 1 x 109.99 110.01 y 0.0 4.0 z -29.9 0
interface 3 face range group 1 x 179.99 180.01 y 0.0 4.0 z -29.9 0
interface 4 face range group 1 x 201.99 202.01 y 0.0 4.0 z -29.9 0
interface 1 prop kn=4e8 ks=4e8 tens=5e3 coh=0.0 fric=20 ;接触面参数
interface 2 prop kn=4e8 ks=4e8 tens=5e3 coh=0.0 fric=20 ;接触面参数
interface 3 prop kn=4e8 ks=4e8 tens=5e3 coh=0.0 fric=20 ;接触面参数
interface 4 prop kn=4e8 ks=4e8 tens=5e3 coh=0.0 fric=20 ;接触面参数
interface 1 maxedge=1
interface 2 maxedge=1
interface 3 maxedge=1
interface 4 maxedge=1
;interface 1 prop kn=4e8 ks=4e8 tens=1e10 sbratio=100
;plot set ba wh
;pl ske interface red blue attach cyan green
;set outp interface_attachment.wmf
;pl ha
;重新定义连续墙参数;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
sel liner id 1 prop isotropic (2.0e10, 0.20) &
cs_nk=4e9 cs_sk=4e9 &
cs_scoh=4e7 cs_scohres=0.0 cs_sfric=0.0 range x 89.9 90.1 y -0.1 4.1 z -30.1 0.1
sel liner id 2 prop isotropic (2.0e10, 0.20) &
cs_nk=4e9 cs_sk=4e9 &
cs_scoh=4e7 cs_scohres=0.0 cs_sfric=0.0 range x 109.9 110.1 y
-0.1 4.1 z -30.1 0.1
sel liner id 3 prop isotropic (2.0e10, 0.20) &
cs_nk=4e9 cs_sk=4e9 &
cs_scoh=4e7 cs_scohres=0.0 cs_sfric=0.0 range x 179.9 180.1 y -0.1 4.1 z -30.1 0.1
sel liner id 4 prop isotropic (2.0e10, 0.20) &
cs_nk=4e9 cs_sk=4e9 &
cs_scoh=4e7 cs_scohres=0.0 cs_sfric=0.0 range x 201.9 202.1 y -0.1 4.1 z -30.1 0.1
;重新定义墙底约束条件;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
def redef_wall_end_link1
node_pnt = nd_head
link_id=100000
loop while node_pnt # null
node_id = nd_id(node_pnt)
xx = nd_pos(node_pnt,2,1)
yy = nd_pos(node_pnt,2,2)
zz = nd_pos(node_pnt,2,3)
link_pnt = nd_link(node_pnt)
dist_x = sqrt((xx-90.0)^2+(zz+30.0)^2)
if dist_x <=dist_tol then
if link_pnt # null then
temp1 = lk_delete(link_pnt)\
link_id = link_id+1
command
sel set link node_tol = dist_tol
sel link id=link_id node_id target zone
sel link attach xdir=rigid ydir=rigid zdir=rigid xrdir=free yrdir=free zrdir=free range id=link_id
endcommand
endif
endif
node_pnt = nd_next(node_pnt)
endloop
end
redef_wall_end_link1
def redef_wall_end_link2
node_pnt = nd_head
link_id=150000
loop while node_pnt # null
node_id = nd_id(node_pnt)
xx = nd_pos(node_pnt,2,1)
yy = nd_pos(node_pnt,2,2)
zz = nd_pos(node_pnt,2,3)
link_pnt = nd_link(node_pnt)
dist_x = sqrt((xx-110.0)^2+(zz+30.0)^2)
dist_tol = 1e-1
if dist_x <=dist_tol then
if link_pnt # null then
if yy < 85.0
temp1 = lk_delete(link_pnt)\
link_id = link_id+1
command
sel set link node_tol = dist_tol
sel link id=link_id node_id target zone
sel link attach xdir=rigid ydir=rigid zdir=rigid xrdir=free yrdir=free zrdir=free range id=link_id
endcommand
endif
endif
endif
node_pnt = nd_next(node_pnt)
endloop
end
redef_wall_end_link2
def redef_wall_end_link3
node_pnt = nd_head
link_id=200000
loop while node_pnt # null
node_id = nd_id(node_pnt)
xx = nd_pos(node_pnt,2,1)
yy = nd_pos(node_pnt,2,2)
zz = nd_pos(node_pnt,2,3)
link_pnt = nd_link(node_pnt)
dist_x = sqrt((xx-180.0)^2+(zz+30.0)^2)
dist_tol = 1e-1
if dist_x <=dist_tol then
if link_pnt # null then
temp1 = lk_delete(link_pnt)\
link_id = link_id+1
command
sel set link node_tol = dist_tol
sel link id=link_id node_id target zone
sel link attach xdir=rigid ydir=rigid zdir=rigid xrdir=free yrdir=free zrdir=free range id=link_id
endcommand
endif
endif
node_pnt = nd_next(node_pnt)
endloop
end
redef_wall_end_link3
def redef_wall_end_link4
node_pnt = nd_head
link_id=250000
loop while node_pnt # null
node_id = nd_id(node_pnt)
xx = nd_pos(node_pnt,2,1)
yy = nd_pos(node_pnt,2,2)
zz = nd_pos(node_pnt,2,3)
link_pnt = nd_link(node_pnt)
dist_x = sqrt((xx-202.0)^2+(zz+30.0)^2)
dist_tol = 1e-1
if dist_x <=dist_tol then
if link_pnt # null then
temp1 = lk_delete(link_pnt)\
link_id = link_id+1
command
sel set link node_tol = dist_tol
sel link id=link_id node_id target zone
sel link attach xdir=rigid ydir=rigid zdir=rigid xrdir=free yrdir=free zrdir=free range id=link_id
endcommand
endif
endif
node_pnt = nd_next(node_pnt)
endloop
end
redef_wall_end_link4
;剑桥模型;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
model cam-clay
;cam-clay模型则不需定义弹性模量(E、G、K)等参数，自动计算
;cam-clay模型中需确定8个模型参数(①-⑧)，手册property中的初始比体积cv(v0)和shear 无须给定
def install_prop
pnt=zone_head
loop while pnt # null
abs_sxx=abs(z_sxx(pnt)) ;|sxx|
abs_syy=abs(z_syy(pnt)) ;|syy|
abs_szz=abs(z_szz(pnt)) ;|szz|
p0=(abs_sxx+abs_syy+abs_szz)/3.0
;cam-clay模型中p、q均须为正值，p0由初应力场确定，故cam-clam定义模型参数前须先已知初应力
p0_effective=p0-z_pp(pnt) ;p0'
;q0=sqrt(((abs_sxx-abs_syy)^2+(abs_syy-abs_szz)^2+(abs_szz-abs_sxx)^2)*0.5)
q0=sqrt(((abs_sxx-abs_syy)^2+(abs_syy-abs_szz)^2+(abs_szz-abs_sxx)^2)*0.5+3.0*((z_sxy(pnt)) ^2+(z_sxz(pnt))^2+(z_syz(pnt))^2))
z_prop(pnt,'mm')=6.0*sin(fai*degrad)/(3.0-sin(fai*degrad)) ;①注三角函数中需将角度转化为弧度
temp1=q0/(z_prop(pnt,'mm')*p0_effective)
pc0=p0_effective*(1.0+temp1^2)*OCR ;先期有效固结压力，用于确定屈服面
v0=1.0+_e0
z_prop(pnt,'cam_cp')=p0_effective ;★重要参数，否则不能正确计算有效应力，提示出错"Mean effective pressure is negative"
z_prop(pnt,'mpc')=pc0 ;②
z_prop(pnt,'poisson')=p_ratio ;③
z_prop(pnt,'lambda')=_lambda ;④
z_prop(pnt,'kappa')=_kappa ;⑤
z_prop(pnt,'mp1')=_mp1 ;⑥
z_prop(pnt,'mv_l')=v0+_lambda*ln(2.0*_cu/(z_prop(pnt,'mm')*_mp1))+(_lambda-_kappa)*l n(2.0) ;⑦
z_prop(pnt,'bulk_bound')=100*40e6 ;⑧
;z_prop(pnt,'bulk_bound')=100*(s_mod+4.0/3.0*s_mod) ;弹性体模上界Kmax
;自动确定Kmax时会出现“property bad”错误提示
;因为弹性上界对计算结果无影响，在不提示Kmax太小的性况下，取值越小计算收敛越快
pnt=z_next(pnt)
endloop
end
set p_ratio=0.25 fai=34.5 _lambda=0.14 _kappa=0.012 _mp1=1e3 _e0=1.2 _cu=10e3 OCR=1.0 ;模型所需参数
install_prop
solve
save model_cam.sav。