- 积分
- 0
- 注册时间
- 2007-9-26
- 仿真币
-
- 最后登录
- 1970-1-1
|
本帖最后由 a02013119 于 2009-12-29 18:59 编辑
为了方便大家FLAC3D的学习,同时也为了simwe这个论坛的浓烈学习气氛,特奉献出自己的原创fish程序。
读书时,发现大家对自己的成果,都有所保留,当然这情由所原,但是这也可能导致学术上的停滞不前,不利于大环境的发展。
基于此,小弟奉献出自己的原创程序,抛砖引玉,希望更多的大牛能够奉献出自己的宝贝,大家为simwe论坛的发展添砖加瓦。营造一个良好的学术气氛。
关于这个程序的使用范围,使用说明中有,原来总结过部分,丢了,现在重新补写了份,更多的功能期待大家去发掘,
关于每个句子的使用,已经添加了注释。
另外,关于这个这个程序,希望大家能够尊重a02013119的劳动果实,不要将之据为己有,去发表文献。
;*******************************************************************************
;********************************2009年5月30日**********************************
;**********************************a02013119编撰*************************************
new ;开始一个新的分析
plot set background white ;背景色设为白色
;plot add block group ;显示所有group
config dyn ;打开动力计算模块
set log on ;建立日志文件,方便查错。
set logfile '日志.txt' ;日志文件名称设为“日志”
;************************************************************
;*************定义建模所需要的参数***************************
def slope_var ;定义边坡形状、结构面产状所需要的参数
ll= 241.44 ;in('('+'输入边坡左侧距离'+'default:'+string(120.0)+'):')
lr= 180.00 ;in('('+'输入边坡右侧距离'+'default:'+string(120.0)+'):')
ld= 120.00 ;in('('+'输入边坡基础高度'+'default:'+string(60.00)+'):')
ls= 100.00 ;in('('+'输入边坡y向距离' +'default:'+string(60.00)+'):')
hd= 0.2000 ;in('('+'输入薄层厚度;无厚度用interface,有厚度用薄层单元'+'default:'+string(0.200)+'):') &
;输数值的时候要注意,为0的话表示采用interface单元,有厚度的话采用薄层单元
dh= 120.00 ;in('('+'输入边坡高度'+'default:'+string(60.0)+'):')
; 输入边坡的倾角参数,一般a_1在45~90之间,如果a_2 = a_3表示单滑面边坡,一般a_2比较小,a_3的数值在a_1左右或大或小皆可。
a_1=45 ;in('('+'输入边坡底部角度'+'default:'+string(15)+'~ '+string(90)+'):')
a_2=30 ;in('('+'输入边坡第二段折线角度'+'<='+string(a_1)+'):')
a_3=50 ;in('('+'输入边坡第一段折线角度'+'在a_1值左右'+'):')
;*******按照高程从低到高 角度依次为a1 a2 a3*********
a1=a_1*degrad ;degrad 是fish 中弧度的意思
a2=a_2*degrad ;为保证滑坡体成型,
a3=a_3*degrad
lt=20.560 ;in('('+'输入滑坡顶端长度'+'<'+string(dh/(1/tan(a2)-1/tan(a1)))+'):')
end
slope_var
;************************************************************
;*************根据输入数据建模并划分网格*********************
def slope_gen ;边坡模型的网格剖分
flag=0 ;定义标记,起始为0
if a2=a3 then ;如果成立,则为单滑面边坡
;*************单滑面边坡*************
xh=lt/(1/tan(a2)-1/tan(a1)) ;单滑面边坡高度是不定的,由此几何关系式推算而来
if xh<=dh then ;单滑面滑坡体成立,必须使得 h < H
ss_x_2_1=ll+lt+dh/tan(a1) ;ss表示 singal slope 单滑面边坡;x表示坐标向;2表示group2;1表示P1,以下雷同。
ss_x_3_1=ll+lt+lr+dh/tan(a1)
ss_x_5_6=ll+xh/tan(a2)
ss_x_6_1=ss_x_5_6-hd*cos(a1)/sin(a1-a2)
ss_x_6_6=ll+hd/sin(a2)
ss_x_7_0=ll+lt
ss_z_4_3=ld+dh
ss_z_5_6=dh+ld-xh
ss_z_6_1=ss_z_5_6+hd*sin(a1)/sin(a1-a2)
flag=1
ssl=30.0 ;边坡剖分的网格尺寸大小
ss_nx_1=round(ll/ssl) ;roud就是四舍五入的意思,由此可得单元的数量
ss_nx_2=round((lt+dh/tan(a1))/ssl)
ss_nx_3=round(lr/ssl)
ss_ny =round(ls/ssl)
ss_nz_1=round(ld/ssl)
ss_nz_4=round(dh/ssl)
ssl7=10.0 ;滑块体剖分的网格尺寸大小
ss_nx_7=round(lt/ssl7)
ss_nz_7=2*ss_nx_2 ;指定薄层单元应用后group 7 z方向的剖分精度
;******************单滑面边坡建模**********************
command
gen zone brick p0(0 0 0) p1(ll 0 0) p2(0 ls 0) p3(0 0 ld) size ss_nx_1 ss_ny ss_nz_1 group 1
gen zone brick p0(ll 0 0) p1(ss_x_2_1 0 0) p2(ll ls 0) p3(ll 0 ld) size ss_nx_2 ss_ny ss_nz_1 group 2
gen zone brick p0(ss_x_2_1 0 0) p1(ss_x_3_1 0 0) p2(ss_x_2_1 ls 0) p3(ss_x_2_1 0 ld) size ss_nx_3 ss_ny ss_nz_1 group 3
gen zone brick p0(0 0 ld) p1(ll 0 ld) p2(0 ls ld) p3(0 0 ss_z_4_3) size ss_nx_1 ss_ny ss_nz_4 group 4
gen zone brick p0(ll 0 ld) p1(ss_x_2_1 0 ld) p2(ll ls ld) p3(ll 0 ss_z_4_3) &
p4(ss_x_2_1 ls ld) p5(ll ls ss_z_4_3) p6(ss_x_5_6 0 ss_z_5_6) &
p7(ss_x_5_6 ls ss_z_5_6) size ss_nx_2 ss_ny ss_nz_4 group 5
;gen zone brick p0(ss_x_2_1 0 ld) p1(ss_x_3_1 0 ld) p2(ss_x_2_1 ls ld) p3(ss_x_5_6 0 ss_z_5_6) &
; p4(ss_x_3_1 ls ld) p5(ss_x_5_6 ls ss_z_5_6) p6(ss_x_3_1 0 ss_z_5_6) &
; p7(ss_x_3_1 ls ss_z_5_6) size ss_nx_3 ss_ny ss_nz_4 group 8 ;开挖第三块
endcommand
if hd =0.0 then ;在单滑面的边坡循环里面,由厚度判断是否用interface单元
command ;薄层厚度为零,故采用interface单元。
gen zone wedge p0(ss_x_7_0 0 ss_z_4_3) p1(ss_x_6_6 0 ss_z_4_3) p2(ss_x_7_0 ls ss_z_4_3) &
p3(ss_x_6_1 0 ss_z_6_1) p4(ss_x_6_6 ls ss_z_4_3)p5(ss_x_6_1 ls ss_z_6_1) &
size ss_nx_7 ss_ny ss_nx_2 group 7
;gen zone brick p0(ss_x_6_1 0 ss_z_6_1) p1(ss_x_3_1 0 ss_z_6_1) p2(ss_x_6_1 ls ss_z_6_1) &
; p3(ss_x_7_0 0 ss_z_4_3) p4(ss_x_3_1 ls ss_z_6_1) p5(ss_x_7_0 ls ss_z_4_3) &
; p6(ss_x_3_1 0 ss_z_4_3) p7(ss_x_3_1 ls ss_z_4_3) size ss_nx_3 ss_ny ss_nx_2 group 10 ;开挖第一块
generate separate 7 ;分离块体7,接触面同一点位置,会产生两个不同的节点
interface 1 wrap 7 5 ;在group 5 和 7 之间设置接触面
;int 1 maxedge 5 ;为了计算精度,可以修改接触面的最大单元尺寸
;generate merge 1e-4 range group 7 any group 10 any ;合并节点
endcommand
ii=out('单滑面边坡,结构面采用interface单元') ;交互输出,人机交流
flag=10
else
command ;如果有厚度,对软弱结构面依然采用薄层单元建模。
gen zone brick p0(ss_x_5_6 0 ss_z_5_6) p1(ss_x_6_1 0 ss_z_6_1) p2(ss_x_5_6 ls ss_z_5_6) &
p3(ll 0 ss_z_4_3) p4(ss_x_6_1 ls ss_z_6_1) p5(ll ls ss_z_4_3) &
p6(ss_x_6_6 0 ss_z_4_3) p7(ss_x_6_6 ls ss_z_4_3) size 1 ss_ny ss_nx_2 group 6
gen zone wedge p0(ss_x_7_0 0 ss_z_4_3) p1(ss_x_6_6 0 ss_z_4_3) p2(ss_x_7_0 ls ss_z_4_3) &
p3(ss_x_6_1 0 ss_z_6_1) p4(ss_x_6_6 ls ss_z_4_3)p5(ss_x_6_1 ls ss_z_6_1) &
size ss_nx_7 ss_ny ss_nz_7 group 7
;对于group7 ,由于设置了薄层单元,可以修改使得它的网格精度是group 5、6的整数倍
attach face plane dip 90.0 dd -30.0 origin ss_x_6_1 0 ss_z_6_1
;当 ss_nz_7=2*ss_nx_2 这一行中出现大于1的整数时才启用上一句,不然就注释掉。
;gen zone brick p0(ss_x_5_6 0 ss_z_5_6) p1(ss_x_3_1 0 ss_z_5_6) p2(ss_x_5_6 ls ss_z_5_6) &
; p3(ss_x_6_1 0 ss_z_6_1) p4(ss_x_3_1 ls ss_z_5_6) p5(ss_x_6_1 ls ss_z_6_1) &
; p6(ss_x_3_1 0 ss_z_6_1) p7(ss_x_3_1 ls ss_z_6_1) size ss_nx_3 ss_ny 1 group 9 ;开挖第二块
;
;gen zone brick p0(ss_x_6_1 0 ss_z_6_1) p1(ss_x_3_1 0 ss_z_6_1) p2(ss_x_6_1 ls ss_z_6_1) &
; p3(ss_x_7_0 0 ss_z_4_3) p4(ss_x_3_1 ls ss_z_6_1) p5(ss_x_7_0 ls ss_z_4_3) &
; p6(ss_x_3_1 0 ss_z_4_3) p7(ss_x_3_1 ls ss_z_4_3) size ss_nx_3 ss_ny ss_nx_2 group 10 ;开挖第一块
endcommand
ii=out('单滑面边坡,结构面采用薄层单元')
flag=11
endif
else
ii=out('单滑面边坡不成立,滑坡高度大于边坡高度')
flag=-1
endif
else
;******************双折线滑面边坡建模**********************
xh=80.00
h3=(xh*(1/tan(a2)-1/tan(a1))-lt)/(1/tan(a2)-1/tan(a3)) ;h3\h2\h1的说明解释请参考边坡图表。
h2=xh-h3
h1=dh-xh ;ds:double slope - 双滑面,其余同单滑面边坡
ds_x_2_1=ll+h3/tan(a3)
ds_x_3_1=ll+h3/tan(a3)+h2/tan(a2)+h1/tan(a1)
ds_x_4_1=ll+lr+h3/tan(a3)+h2/tan(a2)+h1/tan(a1)
ds_x_7_6=ll+h3/tan(a3)+h2/tan(a2)
ds_x_8_6=ll+hd/sin(a3)
ds_x_10_2=ll+lt
ds_x_9_6=ds_x_7_6-hd*cos(a1)/sin(a1-a2)
ds_z_5_3=ld+dh
ds_z_6_6=ld+dh-h3
ds_z_7_6=ld+h1
ds_z_8_1=ld+dh-h3+hd/cos(a3)
ds_z_9_6=ds_z_7_6+hd*sin(a1)/sin(a1-a2)
dsl=20.0 ;网格剖分精度
ds_nx_1=round(ll/dsl)
ds_nx_2=round((h3/tan(a3))/dsl)
ds_nx_3=round((h1/tan(a1)+h2/tan(a2))/dsl)
ds_nx_4=round(lr/dsl)
;用之前询问戴,问是否需要指定剖分精度
ds_nx_6=4*ds_nx_2 ;4表示group 6 剖分精度是 group 2 的4倍
ds_nx_7=2*ds_nx_3 ;2表示group 7 剖分精度是 group 3 的2倍
ds_nz_1=round(ld/dsl)
ds_nz_5=round(dh/dsl)
ds_nx_10=2*ds_nx_7 ;2表示group 10 剖分精度是 group 2 的2倍
ds_nz_10=1*ds_nx_6 ;1表示group 10 剖分精度是 group 6 的1倍
ds_ny =round(ls/dsl)
command
gen zone brick p0(0 0 0) p1(ll 0 0) p2(0 ls 0) p3(0 0 ld) size ds_nx_1 ds_ny ds_nz_1 group 1
gen zone brick p0(ll 0 0) p1(ds_x_2_1 0 0) p2(ll ls 0) p3(ll 0 ld) size ds_nx_2 ds_ny ds_nz_1 group 2
gen zone brick p0(ds_x_2_1 0 0) p1(ds_x_3_1 0 0) p2(ds_x_2_1 ls 0) p3(ds_x_2_1 0 ld) size ds_nx_3 ds_ny ds_nz_1 group 3
gen zone brick p0(ds_x_3_1 0 0) p1(ds_x_4_1 0 0) p2(ds_x_3_1 ls 0) p3(ds_x_3_1 0 ld) size ds_nx_4 ds_ny ds_nz_1 group 4
gen zone brick p0(0 0 ld) p1(ll 0 ld) p2(0 ls ld) p3(0 0 ds_z_5_3) p4(ll ls ld) &
p5(0 ls ds_z_5_3) p6(ll 0 ds_z_5_3) p7(ll ls ds_z_5_3) size ds_nx_1 ds_ny ds_nz_5 group 5
gen zone brick p0(ll 0 ld) p1(ds_x_2_1 0 ld) p2(ll ls ld) p3(ll 0 ds_z_5_3) &
p4(ds_x_2_1 ls ld) p5(ll ls ds_z_5_3) p6(ds_x_2_1 0 ds_z_6_6) &
p7(ds_x_2_1 ls ds_z_6_6) size ds_nx_6 ds_ny ds_nz_5 group 6
gen zone brick p0(ds_x_2_1 0 ld) p1(ds_x_3_1 0 ld) p2(ds_x_2_1 ls ld) &
p3(ds_x_2_1 0 ds_z_6_6) p4(ds_x_3_1 ls ld) p5(ds_x_2_1 ls ds_z_6_6) &
p6(ds_x_7_6 0 ds_z_7_6) p7(ds_x_7_6 ls ds_z_7_6) size ds_nx_7 ds_ny ds_nz_5 group 7
;gen zone brick p0(ds_x_3_1 0 ld) p1(ds_x_4_1 0 ld) p2(ds_x_3_1 ls ld) &
; p3(ds_x_7_6 0 ds_z_7_6) p4(ds_x_4_1 ls ld) p5(ds_x_7_6 ls ds_z_7_6) &
; p6(ds_x_4_1 0 ds_z_7_6) p7(ds_x_4_1 ls ds_z_7_6) size ds_nx_4 ds_ny ds_nz_5 group 11
;如果戴老师指定了group 6 7 是group 2 3 的网格倍数,则下一句是需要的
attach face plane dip 90.0 dd -90.0 origin ds_x_2_1 0 ld range x ll,ds_x_3_1
endcommand
if hd =0.0 then
command ;薄层厚度为零,故采用interface单元。
gen zone brick p0(ll 0 ds_z_5_3) p1(ll ls ds_z_5_3) p2(ds_x_10_2 0 ds_z_5_3) &
p3(ds_x_2_1 0 ds_z_6_6) p4(ds_x_10_2 ls ds_z_5_3) p5(ds_x_7_6 0 ds_z_7_6) &
p6(ds_x_2_1 ls ds_z_6_6) p7(ds_x_7_6 ls ds_z_7_6) size ds_ny ds_nx_10 ds_nz_10 group 10
;gen zone brick p0(ds_x_7_6 0 ds_z_7_6) p1(ds_x_4_1 0 ds_z_7_6) p2(ds_x_7_6 ls ds_z_7_6) &
; p3(ds_x_10_2 0 ds_z_5_3) p4(ds_x_4_1 ls ds_z_7_6) p5(ds_x_10_2 ls ds_z_5_3) &
; p6(ds_x_4_1 0 ds_z_5_3) p7(ds_x_4_1 ls ds_z_5_3) size ds_nx_4 ds_ny ds_nx_2 group 13
generate separate 10
interface 1 wrap 10 6
interface 1 wrap 10 7
;int 1 maxedge 5
;interface 1 prop ks=0.44834308e9 kn=1.210526e9 fric=19.29983041 cohe=0.05e6
;generate merge 1e-4 range group 13 any group 10 any
endcommand
ii=out('双滑面边坡,结构面采用interface单元')
flag=20
else
command
gen zone brick p0(ds_x_2_1 0 ds_z_6_6) p1(ds_x_2_1 0 ds_z_8_1) p2(ds_x_2_1 ls ds_z_6_6) &
p3(ll 0 ds_z_5_3) p4(ds_x_2_1 ls ds_z_8_1) p5(ll ls ds_z_5_3) &
p6(ds_x_8_6 0 ds_z_5_3) p7(ds_x_8_6 ls ds_z_5_3) size 1 ds_ny ds_nx_2 group 8
gen zone brick p0(ds_x_2_1 0 ds_z_6_6) p1(ds_x_7_6 0 ds_z_7_6) p2(ds_x_2_1 ls ds_z_6_6) &
p3(ds_x_2_1 0 ds_z_8_1) p4(ds_x_7_6 ls ds_z_7_6) p5(ds_x_2_1 ls ds_z_8_1)&
p6(ds_x_9_6 0 ds_z_9_6) p7(ds_x_9_6 ls ds_z_9_6) size ds_nx_3 ds_ny 1 group 9
gen zone brick p0(ds_x_2_1 0 ds_z_8_1) p1(ds_x_9_6 0 ds_z_9_6) p2(ds_x_2_1 ls ds_z_8_1) &
p3(ds_x_8_6 0 ds_z_5_3) p4(ds_x_9_6 ls ds_z_9_6)p5(ds_x_8_6 ls ds_z_5_3) &
p6(ds_x_10_2 0 ds_z_5_3) p7(ds_x_10_2 ls ds_z_5_3) size ds_nx_10 ds_ny ds_nz_10 group 10
;gen zone brick p0(ds_x_7_6 0 ds_z_7_6) p1(ds_x_4_1 0 ds_z_7_6) p2(ds_x_7_6 ls ds_z_7_6) &
; p3(ds_x_9_6 0 ds_z_9_6) p4(ds_x_4_1 ls ds_z_7_6) p5(ds_x_9_6 ls ds_z_9_6)&
; p6(ds_x_4_1 0 ds_z_9_6) p7(ds_x_4_1 ls ds_z_9_6) size ds_nx_4 ds_ny 1 group 12
;gen zone brick p0(ds_x_9_6 0 ds_z_9_6) p1(ds_x_4_1 0 ds_z_9_6) p2(ds_x_9_6 ls ds_z_9_6) &
; p3(ds_x_10_2 0 ds_z_5_3) p4(ds_x_4_1 ls ds_z_9_6) p5(ds_x_10_2 ls ds_z_5_3)&
; p6(ds_x_4_1 0 ds_z_5_3) p7(ds_x_4_1 ls ds_z_5_3) size ds_nx_4 ds_ny ds_nz_10 group 13
attach face plane dip 90.0 dd -50.0 origin ds_x_2_1 0 ds_z_8_1 ;如果group 10 网格是左右网格的整数倍时,则启用
endcommand
ii=out('双滑面边坡,结构面采用'+string(hd)+'m薄层单元')
flag=21
endif
endif
end
slope_gen
;****************************************************************
def derive ;根据弹模和泊松比计算K 和 s
s_mod_foundation = y_mod_foundation / (2.0 * (1.0 + p_ratio_foundation))
b_mod_foundation = y_mod_foundation / (3.0 * (1.0 - 2.0 * p_ratio_foundation))
s_mod_up = y_mod_up / (2.0 * (1.0 + p_ratio_up))
b_mod_up = y_mod_up / (3.0 * (1.0 - 2.0 * p_ratio_up))
s_mod_slope = y_mod_slope / (2.0 * (1.0 + p_ratio_slope))
b_mod_slope = y_mod_slope / (3.0 * (1.0 - 2.0 * p_ratio_slope))
s_mod_baoceng = y_mod_baoceng / (2.0 * (1.0 + p_ratio_baoceng))
b_mod_baoceng = y_mod_baoceng / (3.0 * (1.0 - 2.0 * p_ratio_baoceng))
end
set y_mod_foundation = 22e9 p_ratio_foundation = 0.25
set y_mod_up = 12e9 p_ratio_up = 0.25
set y_mod_slope = 12e9 p_ratio_slope = 0.25
set y_mod_baoceng = 1.2e9 p_ratio_baoceng = 0.25
derive
macro foundation 'dens 2850 bulk b_mod_foundation shear s_mod_foundation'
macro up 'dens 2700 bulk b_mod_up shear s_mod_up'
macro slope 'dens 2700 bulk b_mod_slope shear s_mod_slope'
macro baoceng 'dens 2700 bulk b_mod_baoceng shear s_mod_baoceng cohe 0.1e6 fric 28.9 '
macro jiechumian 'ks=0.44834308e9 kn=1.210526e9 fric=37 cohe=0.1e6'
def setup
if flag < 15 then
p_gp=gp_near(ss_x_6_1,0,ss_z_6_1+10)
gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)
;gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp)
;gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp)
gp_disp1 = sqrt(gp_disp)
id1=gp_id(p_gp)
p_gp=gp_near(ss_x_6_1,ls,ss_z_6_1)
gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)
;gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp)
;gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp)
gp_disp2 = sqrt(gp_disp)
id2=gp_id(p_gp)
p_gp=gp_near(ss_x_6_6,0,ss_z_4_3)
gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)
;gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp)
;gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp)
gp_disp3 = sqrt(gp_disp)
id3=gp_id(p_gp)
p_gp=gp_near(ss_x_6_6,ls,ss_z_4_3)
gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)
;gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp)
;gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp)
gp_disp4 = sqrt(gp_disp)
id4=gp_id(p_gp)
else
if flag=20 then
p_gp=gp_near(ds_x_7_6,0,ds_z_7_6)
gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)
;gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp)
;gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp)
gp_disp1 = sqrt(gp_disp)
id1=gp_id(p_gp)
p_gp=gp_near(ds_x_7_6,ls,ds_z_7_6)
gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)
;gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp)
;gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp)
gp_disp2 = sqrt(gp_disp)
id2=gp_id(p_gp)
p_gp=gp_near(ds_x_2_1,0,ds_z_6_6)
gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)
;gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp)
;gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp)
gp_disp1 = sqrt(gp_disp)
id3=gp_id(p_gp)
p_gp=gp_near(ds_x_2_1,ls,ds_z_6_6)
gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)
;gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp)
;gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp)
gp_disp2 = sqrt(gp_disp)
id4=gp_id(p_gp)
p_gp=gp_near(ll,0,ds_z_5_3)
gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)
;gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp)
;gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp)
gp_disp3 = sqrt(gp_disp)
id5=gp_id(p_gp)
p_gp=gp_near(ll,ls,ds_z_5_3)
gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)
;gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp)
;gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp)
gp_disp4 = sqrt(gp_disp)
id6=gp_id(p_gp)
else
p_gp=gp_near(ds_x_9_6,0,ds_z_9_6)
gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)
;gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp)
;gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp)
gp_disp1 = sqrt(gp_disp)
id1=gp_id(p_gp)
p_gp=gp_near(ds_x_9_6,ls,ds_z_9_6)
gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)
;gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp)
;gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp)
gp_disp2 = sqrt(gp_disp)
id2=gp_id(p_gp)
p_gp=gp_near(ds_x_2_1,0,ds_z_8_1)
gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)
;gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp)
;gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp)
gp_disp1 = sqrt(gp_disp)
id3=gp_id(p_gp)
p_gp=gp_near(ds_x_2_1,ls,ds_z_8_1)
gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)
;gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp)
;gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp)
gp_disp2 = sqrt(gp_disp)
id4=gp_id(p_gp)
p_gp=gp_near(ds_x_8_6,0,ds_z_5_3)
gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)
;gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp)
;gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp)
gp_disp3 = sqrt(gp_disp)
id5=gp_id(p_gp)
p_gp=gp_near(ds_x_8_6,ls,ds_z_5_3)
gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)
;gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp)
;gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp)
gp_disp4 = sqrt(gp_disp)
id6=gp_id(p_gp)
endif
endif
end
;setup
def get_gp_maxdisp
gp0_maxdisp = 0
pgp=gp_head
loop while pgp # null
gp_maxdisp = gp_xdisp(pgp)*gp_xdisp(pgp)
gp_maxdisp = gp_maxdisp + gp_ydisp(pgp)*gp_ydisp(pgp)
gp_maxdisp = gp_maxdisp + gp_zdisp(pgp)*gp_zdisp(pgp)
gp_maxdisp = sqrt(gp_maxdisp)
if gp0_maxdisp<gp_maxdisp
gp0_maxdisp=gp_maxdisp
x=gp_xpos(pgp)
y=gp_ypos(pgp)
z=gp_zpos(pgp)
id=gp_id(pgp)
endif
pgp = gp_next(pgp)
endloop
end
;get_gp_maxdisp
def calc_volume
p_z=zone_head
v=0.0
loop while p_z # null
if z_group(p_z)='6'
v=v+z_volume(p_z)
endif
p_z=z_next(p_z)
endloop
end
calc_volume
def get_plast
shearnow = 1
tensionnow = 2
shearpast = 4 ;的state value分别为1、2、4、8
tensionpast= 8 ;具体内容参见fish reference 2.5.3.3
v_shear_now = 0 ;需求的shearnow、tensionnow
v_tension_now = 0 ;shearpast、tensionpast的体积,初始值为0
v_shear_past = 0
v_tension_past = 0 ;flac默认的值
;xx=z_state(find_zone(11779),0) ;通过本步可验证
p_z = zone_head ;单元指针
loop while p_z # null ;开始循环,如果单元指针不为空,执行以下操作
if and(z_state(p_z,0),shearnow) = shearnow then ;判断语句,判断单元如果是shearnow状态
v_shear_now = v_shear_now + z_volume(p_z) ;则把该单元的体积加进v_shear_now
endif
if and(z_state(p_z,0),tensionnow) = tensionnow then ;同上判断语句,判断单元是否为tensionnow状态
v_tension_now = v_tension_now + z_volume(p_z)
endif
if and(z_state(p_z,0),shearpast) = shearpast then ;同上,判断单元是否shearpast
v_shear_past = v_shear_past + z_volume(p_z)
endif
if and(z_state(p_z,0),tensionpast) = tensionpast then
v_tension_past = v_tension_past + z_volume(p_z)
endif
p_z = z_next(p_z) ;单元指针指向下一个单元
endloop
vratio_shear_now=v_shear_now/v ;计算剪切破坏的百分比
vratio_tension_now=v_tension_now/v ;计算张拉破坏的百分比
vratio_shear_past=v_shear_past/v
vratio_tension_past=v_tension_past/v
xii = out('shear_now : ' + string(vratio_shear_now)) ;输出shear_now的体积值
xii = out('tension_now : ' + string(vratio_tension_now)) ;输出tension_now的体积值
xii = out('shear_past : ' + string(vratio_shear_past)) ;输出shear_past的体积值
xii = out('tension_past : ' + string(vratio_tension_past)) ;输出tension_past的体积值
end
;get_plast ;执行函数
;;;运行fish后,执行pri fish就可以显示shearnow(当前剪切破坏),tensionnow(当前拉帐破坏)
;;;shearpast(过去剪切;破坏),tensionpast(过去拉帐破坏)的体积值
;************************************************************************************
def calc_100 ;记录k=1.00时刻的特征点的位移和塑性区破坏比率
xtable(1,1)=0.00
xtable(2,1)=0.00
xtable(3,1)=0.00
xtable(4,1)=0.00
ytable(1,1)=1.00
ytable(2,1)=1.00
ytable(3,1)=1.00
ytable(4,1)=1.00
xtable(5,1)=vratio_shear_now
xtable(6,1)=vratio_tension_now
ytable(5,1)=1.00
ytable(6,1)=1.00
xtable(7,1)=0.00
ytable(7,1)=1.00
end
;calc_100
def ssr ;折减计算
ssx_1=ss_x_3_1-1.0 ;单滑面右侧边界的范围
ssx_2=ss_x_3_1+1.0
dsx_1=ds_x_4_1-1.0 ;双滑面右侧边界的范围
dsx_2=ds_x_4_1+1.0
loop i (1,4) ;折减10次,如若要计算天然状态,只需改成(0,0)
sum=1.00 ;初始值为1.00
n=sum+i*0.05 ;每次折减间隔为0.05,即1.05-1.10-1.15 .....
nn=100.0*n ;
nnn=int(nn) ;数据类型转换
ii=i+1
name1 = 'This is run number ' + string(n)
name2 = 'flac' + string(nnn) + '.sav'
coh1=0.1e6/n ;粘聚力折减
fri1=(atan((tan(37*degrad))/n))*180/pi ;摩擦系数折减
if flag > 5 then
if flag < 20 then
if flag=10 then
command
model null ;清空上一次计算的保留场
model elastic
prop foundation range group 1
prop foundation range group 2
prop foundation range group 3
prop up range group 4
prop up range group 5
prop slope range group 7
;prop dens 2700 bulk 2.314815e8 shear 0.94697e8 range group 8
;prop dens 2700 bulk 2.314815e8 shear 0.94697e8 range group 10
interface 1 jiechumian
set grav 0 0 -9.81
fix x y z range z -0.1 0.1
fix x range x -1 1
fix x range x ssx_1 ssx_2 ;定义边界范围
fix y
his unbal
plot show
plot set background white
plot add contour zdisp outline on
set dyn off ;关闭动力计算模块
set echo off
solve
save flac100
calc_100 ;调用clac_100函数,计算flac100的位移值
title name1
ini xvel=0 yvel=0 zvel=0 xdisp=0 ydisp=0 zdisp=0
interface 1 prop ks=0.44834308e9 kn=1.210526e9 fric=fri1 cohe=coh1
solve
save name2
endcommand
;get_gp_maxdisp
setup
xtable(1,ii)=gp_disp1 ;把这个折减系数下的特征点1的位移赋值给table 1 中x列的第ii行
xtable(2,ii)=gp_disp2
xtable(3,ii)=gp_disp3
xtable(4,ii)=gp_disp4
ytable(1,ii)=n
ytable(2,ii)=n
ytable(3,ii)=n
ytable(4,ii)=n
xtable(7,ii)=gp_maxdisp
ytable(7,ii)=n
else
command
model null
model elastic
prop foundation range group 1
prop foundation range group 2
prop foundation range group 3
prop up range group 4
prop up range group 5
prop slope range group 7
;prop dens 2700 bulk 2.314815e8 shear 0.94697e8 range group 8
;prop dens 2700 bulk 2.314815e8 shear 0.94697e8 range group 9
;prop dens 2700 bulk 2.314815e8 shear 0.94697e8 range group 10
model mohr range group 6 any
prop baoceng range group 6 any
set grav 0 0 -9.81
fix x y z range z -1 1
fix x range x -1 1
fix x range x ssx_1 ssx_2
fix y
hist unbal
plot show
plot add hist 1
plot set background white
plot add contour zdisp outline on
set dyn off
set echo off
solve
save flac100
get_plast
calc_100 ;调用clac_100函数,计算flac100的位移值
title name1
ini xvel=0 yvel=0 zvel=0 xdisp=0 ydisp=0 zdisp=0
prop dens 2700 bulk 0.730111e9 shear 0.24337e9 cohesion coh1 friction fri1 range group 6
solve
save name2
endcommand
setup
get_plast
;get_gp_maxdisp
xtable(1,ii)=gp_disp1
xtable(2,ii)=gp_disp2
xtable(3,ii)=gp_disp3
xtable(4,ii)=gp_disp4
xtable(5,ii)=vratio_shear_now
xtable(6,ii)=vratio_tension_now
xtable(7,ii)=gp0_maxdisp
ytable(1,ii)=n
ytable(2,ii)=n
ytable(3,ii)=n
ytable(4,ii)=n
ytable(5,ii)=n
ytable(6,ii)=n
ytable(7,ii)=n
endif
else
if flag=20 then
command
model null
model elastic
prop foundation range group 1
prop foundation range group 2
prop foundation range group 3
prop foundation range group 4
prop up range group 5
prop up range group 6
prop up range group 7
prop baoceng range group 10
;prop dens 2700 bulk 2.314815e8 shear 0.94697e8 range group 11
;prop dens 2700 bulk 2.314815e8 shear 0.94697e8 range group 13
interface 1 jiechumian
set grav 0 0 -9.81
fix x y z range z -1 1
fix x range x -1 1
fix x range x dsx_1 dsx_2
fix y
plot show
plot set background white
plot add contour zdisp outline on
set dyn off
set echo off
solve
save flac100
calc_100
title name1
ini xvel=0 yvel=0 zvel=0 xdisp=0 ydisp=0 zdisp=0
interface 1 prop ks=0.44834308e9 kn=1.210526e9 fric=fri1 cohe=coh1
solve
save name2
endcommand
setup
;get_gp_maxdisp
xtable(1,ii)=gp_disp1
xtable(2,ii)=gp_disp2
xtable(3,ii)=gp_disp3
xtable(4,ii)=gp_disp4
xtable(7,ii)=gp_maxdisp
ytable(1,ii)=n
ytable(2,ii)=n
ytable(3,ii)=n
ytable(4,ii)=n
ytable(7,ii)=n
else
command
model elastic
group 89 range group 8 any group 9 any
prop foundation range group 1
prop foundation range group 2
prop foundation range group 3
prop foundation range group 4
prop up range group 5
prop up range group 6
prop up range group 7
prop slope range group 10
;prop dens 2700 bulk 2.314815e8 shear 0.94697e8 range group 11
;prop dens 2700 bulk 2.314815e8 shear 0.94697e8 range group 12
;prop dens 2700 bulk 2.314815e8 shear 0.94697e8 range group 13
model mohr range group 89 any
prop baoceng range group 89 any
set grav 0 0 -9.81
fix x y z range z -1 1
fix x range x -1 1
fix x range x dsx_1 dsx_2
fix y
plot show
plot set background white
plot add contour zdisp outline on
set dyn off
set echo off
solve
save flac100
get_plast
calc_100
title name1
ini xvel=0 yvel=0 zvel=0 xdisp=0 ydisp=0 zdisp=0
prop dens 2700 bulk 0.730111e9 shear 0.24337e9 cohesion coh1 friction fri1 range group 89 any
solve
save name2
endcommand
setup
get_plast
xtable(1,ii)=gp_disp1
xtable(2,ii)=gp_disp2
xtable(3,ii)=gp_disp3
xtable(4,ii)=gp_disp4
xtable(5,ii)=vratio_shear_now
xtable(6,ii)=vratio_tension_now
xtable(7,ii)=gp_maxdisp
ytable(1,ii)=n
ytable(2,ii)=n
ytable(3,ii)=n
ytable(4,ii)=n
ytable(5,ii)=n
ytable(6,ii)=n
ytable(7,ii)=n
endif
endif
endif
endloop
end
ssr
def disp ;给位移折减曲线图形添加标示
array ve(3)
ve(1) = -200.00 ;x向坐标
ve(2) = 0.00 ;y向坐标
ve(3) = 240.0 ;z向坐标;三向定位
oo = set_fontsize(0.5) ;设置字体大小
oo = draw_string(ve,'折减系数-K')
ve(1) = 300.0
ve(3) = -180.00
oo = set_fontsize(0.5)
oo = draw_string(ve,'x-disp');
end ;给图形添加标示
def state ;给塑性区折减曲线图形添加标示
array vec(3)
vec(1) = -200.00
vec(2) = 0.00
vec(3) = 240.0
oo = set_fontsize(0.5) ;设置字体大小
oo = draw_string(vec,'折减系数-K')
vec(1) = 300.0
vec(3) = -180.00
oo = set_fontsize(0.5)
oo = draw_string(vec,'塑性区百分比');
end ;给图形添加标示
table 1 name '特征点1位移与折减系数关系曲线' ;定义table的名称
table 2 name '特征点2位移与折减系数关系曲线'
table 3 name '特征点3位移与折减系数关系曲线'
table 4 name '特征点4位移与折减系数关系曲线'
table 5 name '剪切破坏区与折减系数关系曲线'
table 6 name '张拉破坏区与折减系数关系曲线'
table 7 name '最大位移与折减系数关系曲线'
plot creat '位移与折减系数关系曲线’ ;创建一个新图形窗口
plot set background white
plot add table 1 2 3 4 both
plo add fish disp black
plot creat '最大位移与折减系数关系曲线’
plot set background white
plot add table 7 line
;开始进行动力计算
set dyn on
set large
ini xvel=0 yvel=0 zvel=0 xdisp=0 ydisp=0 zdisp=0
ini state = 0
free x y z
;采用刚性地基输入法**************************************
table 1 read vel.txt
range name bottom z=-.1 .1
apply xvel=0.859852276 hist table 1 range bottom
apply yvel=0.51054291 hist table 1 range bottom
fix z ran z -1 1
set dyn damp local 0.15
;*********************************************************
;采用柔性地基输入法**************************************
;range name bottom z=-.1 .1
;apply nquiet squiet dquiet range plane norm 0,0,1 bottom
;table 1 read vel.txt
;apply dstress=6022488 hist table 1 range bottom
;apply sstress=-3575892 hist table 1 range bottom
;apply nvel 0 plane norm 0,0,1 range bottom
;set dyn damp local 0.15
;*********************************************************
apply ff
set dyn time =0
hist reset
hist dytime
;添加动力的检测点
plot show
plot set rotation 20 0 270
plot set mag 1.25
plot add contour xdisp outline on
hist n=100
set dyn multi on
def su
old_time=clock
end
su
def tim
tim=0.01*(clock-old_time)
end
solve age 30.000
save 30
pri tim
set log off |
评分
-
1
查看全部评分
-
|