找回密码
 注册
Simdroid-非首页
查看: 9881|回复: 387

[命令/FISH] 毕业后倾情奉献边坡建模、计算、后处理综合fish程序(原创

[复制链接]
发表于 2009-12-8 19:55:46 | 显示全部楼层 |阅读模式 来自 江苏南京
本帖最后由 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

查看全部评分

 楼主| 发表于 2009-12-29 18:50:56 | 显示全部楼层 来自 江苏南京
Simdroid开发平台
;*******************************************************************************
;********************************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
回复 9 不支持 0

使用道具 举报

 楼主| 发表于 2009-12-9 08:41:53 | 显示全部楼层 来自 江苏南京
本帖最后由 a02013119 于 2009-12-9 08:43 编辑

该程序将某典型边坡建模(尺寸)、剖分(单元尺寸精度)、计算、折减计算、后处理画图,综合为一体,只要输入计算参数,最后直接得出折减安全系数。很方便。

看人气,后有更多奉献。
回复 2 不支持 1

使用道具 举报

发表于 2009-12-17 21:50:35 | 显示全部楼层 来自 江苏南京
送人玫瑰手留余香,手留余香心情快乐啊,愿a02013119天天快乐!!

点评

楼主真无私,顶  发表于 2013-5-19 09:37
楼主真的很给力啊  发表于 2013-5-13 19:32
回复 1 不支持 0

使用道具 举报

发表于 2009-12-8 20:35:40 | 显示全部楼层 来自 湖南长沙
谢谢无私奉献,学一下
回复 1 不支持 0

使用道具 举报

发表于 2009-12-8 21:24:30 | 显示全部楼层 来自 江苏徐州
楼主真是好人啊!!!
回复 不支持

使用道具 举报

发表于 2009-12-8 23:13:55 | 显示全部楼层 来自 江苏南京
谢谢楼主分享。
回复 不支持

使用道具 举报

发表于 2009-12-8 23:16:52 | 显示全部楼层 来自 重庆大学
论坛高手真多,建议版本主加分啊
回复 不支持

使用道具 举报

发表于 2009-12-9 10:01:29 | 显示全部楼层 来自 北京
像楼主致敬,真希望哪天也能将这个语言应用得挥洒自如
回复 不支持

使用道具 举报

发表于 2009-12-9 10:52:47 | 显示全部楼层 来自 山东烟台
像楼主致敬,真希望哪天也能将这个语言应用得挥洒自如
回复 不支持

使用道具 举报

发表于 2009-12-9 11:35:28 | 显示全部楼层 来自 江苏徐州
谢谢楼主
像您这样的为了大家更好学习进步的人值得大家的敬佩
回复 不支持

使用道具 举报

发表于 2009-12-9 13:18:31 | 显示全部楼层 来自 安徽淮南
精神可嘉 呵呵 支持原创
回复 不支持

使用道具 举报

发表于 2009-12-9 18:20:03 | 显示全部楼层 来自 四川大学
谢谢楼主,正在自学,有很多东西不懂,也很难找到资料学习,难啊~~~~~~
回复 不支持

使用道具 举报

发表于 2009-12-9 18:46:42 | 显示全部楼层 来自 湖北武汉
支持原创性。。。。
回复 不支持

使用道具 举报

发表于 2009-12-9 22:13:59 | 显示全部楼层 来自 四川成都
谢谢楼主的无私奉献,学习ing
回复 不支持

使用道具 举报

发表于 2009-12-9 23:08:47 | 显示全部楼层 来自 湖北武汉
nice ok
不错的哦 可以好好的学习一下
回复 不支持

使用道具 举报

发表于 2009-12-10 09:49:17 | 显示全部楼层 来自 广东深圳
新手入门,感谢分享
回复 不支持

使用道具 举报

发表于 2009-12-10 16:00:39 | 显示全部楼层 来自 广东江门
谢谢
楼主!!!
回复 不支持

使用道具 举报

发表于 2009-12-14 01:35:55 | 显示全部楼层 来自 上海交通大学
感谢楼主分享
回复 不支持

使用道具 举报

发表于 2009-12-14 11:19:09 | 显示全部楼层 来自 四川成都
楼主在贴两张图上来看看就完美了!十分感谢楼主的无私奉献!!!
回复 不支持

使用道具 举报

发表于 2009-12-14 13:31:02 | 显示全部楼层 来自 甘肃兰州
感谢楼主的无私奉献!!!
回复 不支持

使用道具 举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

Archiver|小黑屋|联系我们|仿真互动网 ( 京ICP备15048925号-7 )

GMT+8, 2024-3-29 02:31 , Processed in 0.070171 second(s), 15 queries , Gzip On, MemCache On.

Powered by Discuz! X3.5 Licensed

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表