自编强度折减法算边坡降雨入渗的稳定性，速度矢量是错的，求高手指点

tanghuiyun

new
;建立网格模型
gen zone brick  &
p0 0 0 0 p1 2 0 0 p2 0 0.5 0 p3 0 0 3 size 3 1 3
gen zone brick  &
p0 2 0 0 p1 20 0 0 p2 2 0.5 0 p3 2 0 3 &
size 17 1 3 ratio 1.03 1 1
gen zone brick  &
p0 2 0 3 p1 20 0 3 p2 2 0.5 3 p3 12 0 13 &
p4 20 0.5 3 p5 12 0.5 13 p6 20 0 13 &
p7 20 0.5 13 size 17 1 17 ratio 1.03 1 1
;====================================
def SSR
;定义有关参数及循环终止条件
ait1=0.02
k11=1.0
k12=2.0
ks=(k11+k12)/2
loop while (k12-k11)>ait1
  coh1=7.5e4/ks
  fri1=atan(tan(18*pi/180)/ks)*180/pi
  dila1=18.0
  ten1=-5.5e5
  K1=5.1e10
  G1= 2.35e10
;折减实现过程
  command
    model null
    config fluid
    set fluid off
    pl bl g
    model elas
    prop bulk 5.1e12 shear 2.35e12
    model fl_iso
    prop poros 0.5 perm 1.34e-11
    ini fden 1000
    ini fmod 2e9 ftens -1e-3
    apply nstress -3e4 grad 0 0 10e3 range x -0.1 0.1 y 0 3
;设置边界条件
fix x y z range z -0.1 0.1
fix x range x 19.9 20.1
fix x range x -0.1 0.1
fix y
ini dens 2000
set gravity 0 0 -10.0
history unb
solve
save lizi2.sav
;位移和速度归零
ini xdisp 0
ini ydisp 0
ini zdisp 0
ini xvel 0
ini yvel 0
ini zvel 0
;...fluid flow model...
    model fl_iso
    prop poros 0.5  
    ini fden 1000
    ini fmod 6.7e3 ftens -1e-3
    water density 1000
    table 1 name hanshuilv
    table 2 name jizhixili
    table 3 name shentouxishu
    def ana_sol
      n_max=20
      aerfa=0.006
      c_a=1/aerfa
      c_n=1.45
      c_m=1-1/c_n
      sita_r=0.15
      sita_s=0.5
      perm_s=1e-11
      pnt=gp_head
      loop while pnt # null
         sita0=0.15
         sita=0.0
         n=0     
         loop while n < n_max
            n=n+1
            sita=sita0+n*0.035
            if sita < sita_s
              command
                apply discharge 2.23e-6 range x 2 20 y 0 0.5 z 3 13
              endcommand
            else
              sita=sita_s
              command
                apply nstress -5e2 range x 2 20 y 0 0.5 z 3 13
              endcommand
            endif
            sita1=sita-sita_r
            sita2=sita_s-sita_r
            h_xili=c_a*[(sita1/sita2)^(-1/c_m)-1]^(1-c_m)   
            a_h=aerfa*h_xili
            an_h=a_h^(c_n-1)
            a1_h=1+(a_h)^c_n
            k_perm=perm_s*[1-an_h*a1_h^(-c_m)]^2/a1_h^(c_m/2)
            command
              prop perm k_perm
            endcommand
            table(1,n)=sita
            table(2,n)=h_xili
            table(3,n)=k_perm
         endloop
        pnt=gp_next(pnt)      
      endloop  
    end
    m m
    prop bulk K1 shear G1 coh coh1 fric fri1 tens ten1 dila dila1
    apply nstress -3e4 grad 0 0 10e3 range x -0.1 0.1
;设置边界条件
    fix x y z range z -0.1 0.1
    fix x range x 19.9 20.1
    fix x range x -0.1 0.1
    fix y
;初始化材料密度
    ini dens 2000
   
    set fluid on
    set mech ratio 9.8e-6
    ana_sol
    solve age 43200
    save lizi2_12.sav
    set fluid off
  endcommand
;二分法实现过程
   if mech_ratio<1.0e-5
      k11=ks
      k12=k12
   else
      k12=ks
      k11=k11
   endif
     ks=(k11+k12)/2
endloop
;计算结果的保存
fosfile0='_fos12'+'.sav'
command
  save fosfile0
endcommand
end
SSR
pr ks
[img]file:///D:/Program%20Files/Tencent/QQ/Users/282173282/Image/O6Q(C[`S8ZHX]QYVRJ9ICWQ.jpg[/img] 速度矢量往上走了，孔压也不对，怎么会这样，高手请指教

tanghuiyun

有没有人知道啊？急啊？

fangweino1

运行出错啊
死循环

eagleming

连个图都没有的！！！！

小刚1987

正在学，不知道

YHPXUEXI

不知道楼主解决了吗，正在学习