强度折减法计算边坡稳定性，位移云图数值非常大

lchyn2006

初学者求助，各位大神，我用FLAC3D进行计算边坡稳定性，最终求得边坡安全系数为1.12，是稳定的。
但是通过 plot add cont disp命令查看位移云图却发现位移非常大，最大的位移已经达到几十米远，困扰我很久了。
不知道是我的参数设置不对，还是命令流有问题，请各位大神帮忙！感激不尽！
计算方法：强度折减法，二分法
边坡参数：长290m，宽260m，高130m。
岩石指标：coh=12380pa,fric=20°,tens=1.0e6pa, shear=3e7pa,
bulk=1e8pa, den=0.002e6t/m3.
以下是命令流
*********************************************************
n
set log on
set logfile ceshi.log
restore beikuangdamoxing.sav


;初始应力场生成
m elas
def chushi
jian=yang/(2*(1+passion))
tiji=yang/(3*(1-2*passion))
end
;
;定义材料属性
def cailiao1
command
prop den 0.002e6 shear 3e10 bulk 1e11 ran group 1
endcommand
end
cailiao1


;设置边界条件
fix x y z range Z -0.1 0.1
fix x ;range x -0.1 1.3
;fix x range x 248.5 249.8
fix y range y -0.1 0.1
fix y range y 290.2 290.3

set grav 0  0  -10.0
solve
ini xdisp 0 ydisp 0 zdisp 0
ini xvel 0 yvel 0 zvel 0

;单元自重应力写入文件'graStress.dat'
def graStressInfo
  array buf(6)
  pnt=zone_head
  loop while pnt # null
    ;应力顺序sxx,syy,szz,sxy,sxz,syz
    buf(1)=string(z_sxx(pnt))
    buf(2)=string(z_syy(pnt))
    buf(3)=string(z_szz(pnt))
    buf(4)=string(z_sxy(pnt))
    buf(5)=string(z_sxz(pnt))
    buf(6)=string(z_syz(pnt))
    status = write(buf,6)
    pnt=z_next(pnt)
  end_loop
end
def wrGraStress
  status = close
  status = open('graStress.dat',1,1)
  if status = 0 then
    graStressInfo
    status = close
  end_if
end
wrGraStress


;定义调用应力状态的数组，以免重复定义
def defArray
  array buu(6)
end

;********************************************** ;自定义强度折减法
def SSR  
;=====================================
;定义有关参数及循环终止条件
ait1=0.02      
k11=1.0        
k12=2.0      
ks=(k11+k12)/2  

loop while (k12-k11)>ait1

coh11=12380/ks;slope
fri11=(atan((tan(20*pi/180.0))/ks))*180.0/pi



;=====================================
;折减的实现过程   
command
model null
;用来调用单元自重应力文件'graStress.dat'
defArray
def graStressInfo
  pnt=zone_head
  loop while pnt # null
    ;应力顺序sxx,syy,szz,sxy,sxz,syz
    status = read(buu,6)
    z_sxx(pnt)=float(buu(6))
    z_syy(pnt)=float(buu(1))
    z_szz(pnt)=float(buu(5))
    z_sxy(pnt)=float(buu(2))
    z_sxz(pnt)=float(buu(3))
    z_syz(pnt)=float(buu(4))
    pnt=z_next(pnt)
  end_loop
end
def reGraStress
  status = close
  status = open('graStress.dat',0,1)
  if status = 0 then
    graStressInfo
    status = close
  end_if
end
ini xdisp 0 ydisp 0 zdisp 0
ini xvel 0 yvel 0 zvel 0

;塑性阶段求解
model mohr   

def cailiao11
command
prop shear 3e7 bulk 1e8 ran group 1
prop coh=coh11 fric=fri11 tens=1.0e6 ran group 1
ini den 0.002e6 ran group 1

endcommand
end
cailiao11

;设置边界条件
fix x y z range Z -0.1 0.1
fix x ;range x -0.1 1.3
;fix x range x 248.5 249.8
fix y range y -0.1 0.1
fix y range y 290.2 290.3
set grav 0  0  -10.0

print ks
set mech ratio 9.8e-6     
solve step 6000   
endcommand  

;二分法的实现过程   
if mech_ratio<1.0e-5     
k11=ks     
k12=k12   
else      
k12=ks     
k11=k11   
endif

ks=(k11+k12)/2
endloop  
;=====================================
;计算结果的保存
fosfile0='_fos'+'.sav'
command   
save fosfile0
endcommand
end  
;**********************************************
;程序执行及结果显示
SSR  
pr ks
set log off

chenpf09

你的边界条件那为什么有  ;   呀？？？
强度折减应该只是为了获得稳定系数的吧。

lchyn2006

哦，”；“是我调整边界条件看看有什么影响的。chenpf09您说的对，主要是为了获得稳定系数，可能是我的体积模量和剪切模量取的有问题吧，谢谢您的解答！

走起GOGO

好东西啊

zhaoyigenv

想问下读入命令确定能把应力读进去吗？