模拟岩块单轴加载，为何出现拉应力？

hanxianmin

模拟一个1×1×1m的岩块在无侧压条件下的单轴抗压强度，顶部采用小位移速率加载，发现岩块顶部区域的syy为拉应力（即，正值），并且随着时间步计算应力一直没变化，经过调整仍无法避免这一和经验相矛盾的现象，实在想不出是什么原因，恳请各位大侠予以解答、指正。
源代码如下：
new
set plot pcx
set pcx on
set autoname on
set pcx filename 1.pcx
plot open
movie step on 200
;建模
poly brick 0,1 0,1 0,1
plot hold ax
;取截面观测竖直应力分布
plot xsec dip 90 dd 180 syy ax
plot zone on
;剖分网格
gen edge 0.1
;设置材料力学模型
change cons 2
change jcons 1

;材料属性
def cal_mod
  s_mod = y_mod / (2.0 * (1.0 + p_ratio))
  b_mod = y_mod / (3.0 * (1.0 - 2.0 * p_ratio))
end
set y_mod = 7e10 p_ratio = 0.31
cal_mod

change mat 1
change jmat 1
prop mat=1 dens 3052 k=b_mod g=s_mod bfr=43.5 bc=7e6 bt=9e6

;边界条件
fix 0,1 -0.1,0.1 0,1
bound 0,1 0.99,1.01 0,1 yvel=-1e-5
bound 0,1 0,1 -0.1,0.1 zvel=0.0
bound 0,1 0,1 0.9,1.1 zvel=0.0
;设置重力加速度

grav 0.0 -9.8 0.0
;计算顶面节点平均应力
def cal_str
sum=0.0
  bi=block_head
  m_con=0
  loop while bi # 0
   zi=b_zone(bi)
     loop while zi # 0     
       if z_y(zi)>0.97 then     
         sum=sum+z_syy(zi)
         m_con=m_con+1     
       endif
       zi=z_next(zi)
     end_loop
  bi=b_next(bi)
end_loop
  if m_con # 0 then  
    cal_str=sum/m_con
  end_if
end
;
def cal_strain
  gi=gp_near(0.5,1.0,0.5)
  cal_strain=gp_ydis(gi)
end
;平均应力和位移变化过程控制
hist n 1
hist cal_strain
hist cal_str
;计算
step 2000

plot hist 2 vs 1 xr  hold
;保存文件
save s9.sav
ret

iggcasmcf

请教：既然是单轴抗压，为什么还要把z方向的速度设置为0？