应力的计算

jinbo1122

; FNAME: Shear_box.dat

; the units used in this simulation: mass(kg), length(m), time

(step),density(kg/m^3)

; stiffness(N/m), force(N), pressure(N/m^2)
new
set random
call fishcall.fis
  macro pt1 '6.0 0.0 0.0'
  macro pt2 '6.0 6.0 0.0'
  macro pt3 '0.0 6.0 0.0'
  macro pt4 '0.0 0.0 0.0'
  macro pt5 '6.0 0.0 3.0'
  macro pt6 '6.0 6.0 3.0'
  macro pt7 '0.0 6.0 3.0'
  macro pt8 '0.0 0.0 3.0'
  macro pt9 '6.0 0.0 3.01'
  macro pt10 '6.0 6.0 3.01'
  macro pt11 '0.0 6.0 3.01'
  macro pt12 '0.0 0.0 3.01'
  macro pt13 '6.0 0.0 6.01'
  macro pt14 '6.0 6.0 6.01'
  macro pt15 '0.0 6.0 6.01'
  macro pt16 '0.0 0.0 6.01'
  macro pt17 '6.0 -3.0 3.0'
  macro pt18 '0.0 -3.0 3.0'
  macro pt19 '6.0 9.0 3.0'
  macro pt20 '0.0 9.0 3.0'
  macro pt21 '5.95 0.05 6.01'
  macro pt22 '0.05 0.05 6.01'
  macro pt23 '0.05 5.95 6.01'
  macro pt24 '5.95 5.95 6.01'
  macro wall_s_stiff '1e16'
  macro wall_n_stiff '1e16'
  macro wall_fric  '0.3'
  macro wall_props 'kn wall_n_stiff ks wall_s_stiff fric wall_fric'
  wall id=1 wall_props face pt1 pt2 pt3 pt4
  wall id=2 wall_props face pt2 pt6 pt7 pt3  
  wall id=3 wall_props face pt3 pt7 pt8 pt4
  wall id=4 wall_props face pt4 pt8 pt5 pt1
  wall id=5 wall_props face pt1 pt5 pt6 pt2
  wall id=6 wall_props face pt6 pt19 pt20 pt7
  wall id=7 wall_props face pt5 pt8 pt18 pt17
  wall id=8 wall_props face pt9 pt13 pt14 pt10
  wall id=9 wall_props face pt14 pt15 pt11 pt10  
  wall id=10 wall_props face pt15 pt16 pt12 pt11
  wall id=11 wall_props face pt13 pt9 pt12 pt16
  wall id=12 wall_props face pt13 pt16 pt15 pt14
  delete wall 12
  wall id=13 wall_props face pt21 pt22 pt23 pt24
  


  plot add wall red id on
  plot set rotation (30, 0, 60)
  plot add ball yellow
  plot show

  gen id=1,5000 rad 0.03 0.05 x 0 6 y 0 6 z 0 6.01
  property density 2000 kn 1e7 ks 1e7
  ini rad mul 4.5
  set dt dscale
  set grav 0 0 -9.81
  cycle 100  plot add wall red id on
  plot set rotation (30, 0, 60)
  plot add ball yellow
  plot show
def wall_addr
wadd1 = find_wall(1)
wadd2= find_wall(2)
wadd3 = find_wall(3)
wadd4 = find_wall(4)
wadd5 = find_wall(5)
wadd6 = find_wall(6)
wadd7 = find_wall(7)
wadd8 = find_wall(8)
wadd9= find_wall(9)
wadd10= find_wall(10)
wadd11= find_wall(11)
wadd12= find_wall(12)
wadd13= find_wall(13)
end
;-----------------------------------
def get_ss
wadd1 = find_wall(1)
wadd2= find_wall(2)
wadd3 = find_wall(3)
wadd4 = find_wall(4)
wadd5 = find_wall(5)
wadd6 = find_wall(6)
wadd7 = find_wall(7)
wadd8 = find_wall(8)
wadd9= find_wall(9)
wadd10= find_wall(10)
wadd11= find_wall(11)
wadd12= find_wall(12)
wadd13= find_wall(13)
ydisp1=w_x(wadd4)
ydisp2=w_y(wadd3)
xdif = w_x(wadd8) - w_x(wadd10)
ydif = w_y(wadd11)
zdif = w_z(wadd13)-w_z(wadd1)
wyforce2=w_yfob(wadd2)
wyforce8=w_yfob(wadd8)
wyforce10=w_yfob(wadd10)
wyforce11=w_yfob(wadd11)
wyforce9=w_yfob(wadd9)
wxforce2=w_xfob(wadd2)
wxforce8=w_xfob(wadd8)
wxforce10=w_xfob(wadd10)
wxforce11=w_xfob(wadd11)
wxforce9=w_xfob(wadd9)
shearforce=0.5*(w_yfob(wadd2)-w_yfob(wadd11))
new_xwidth = 6.0 + xdif
new_ywidth = 6.0 + ydif
new_height = 0.5*6.0 + zdif
wsxx = 0.5*(w_xfob(wadd8) - w_xfob(wadd10)) / (new_ywidth * new_height)
wsyy = w_yfob(wadd11) / (new_xwidth * new_height)
wszz = w_zfob(wadd13)/ (new_xwidth * new_ywidth)
wexx = 2.0 * xdif / (6.0 + new_xwidth)
weyy = 2.0 * ydif / (6.0 + new_ywidth)
wezz = 2.0 * zdif / (0.5*6.0 + new_height)
wevol = wexx + weyy + wezz
end
;-----------------------------------------------------------
get_ss
print wadd4
print w_x(wadd4)
;end-------------------------------------------




这是我的程序，理论上没有什么错误，但是w_x( )之中有问题，请高手指导