求帮助

jinbo1122

这是我的一个模拟直剪实验的程序，怎么都调试不通，希望高人指点。
; FNAME: Shear_box.dat

; the units used in this simulation: mass(kg), length(m), time(step),density(kg/m^3) , sxxreq=100000(pa);

; stiffness(N/m), force(N), pressure(N/m^2)

set random

  new

  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

  gen id=1,5000 rad 0.03 0.050 x 0 6 y 0 6 z 0 6.01

  property density 2000 kn 1e8 ks 1e8

  plot add wall red id on

  plot set rotation (30, 0, 60)

  plot add ball yellow

  plot show

  ini rad mul 4.5

  set dt dscale

  cycle 100

  measure id 1 x 0.25 y 0.25 z 0.25 rad 0.24

  print measure 1

  plot add wall red id on

  plot set rotation (30, 0, 60)

  plot add ball yellow

  plot show
     
; cycle 100


   
def wall_addr

w1 = find_wall(1)

w2 = find_wall(2)

w3 = find_wall(3)

w4 = find_wall(4)

w5 = find_wall(5)

w7 = find_wall(7)

w8 = find_wall(8)

w9 = find_wall(9)

w10 = find_wall(10)

w12= find_wall(12)

w13= find_wall(13)

end

;---------------------------------------------
def make_ss ;

ydisp1=w_y(w4)

ydisp2=w_y(w3)

xdif = w_x(w8) - w_x(w10)

ydif = w_y(w11)

zdif = w_z(w13)

wyforce2=w_yfob(w2)

wyforce8=w_yfob(w8)

wyforce10=w_yfob(w10)

wyforce11=w_yfob(w11)

wyforce9=w_yfob(w9)

wxforce2=w_xfob(w2)

wxforce8=w_xfob(w8)

wxforce10=w_xfob(w10)

wxforce11=w_xfob(w11)

wxforce9=w_xfob(w9)

shearforce=0.5*(wyforce2-wyforce11)

new_xwidth = 6.0 + xdif

new_ywidth = 6.0 + ydif

new_height = 0.5*6.0 + zdif

wsxx = 0.5*(w_xfob(w8) - w_xfob(w10)) / (new_ywidth * new_height)

wsyy = w_yfob(w11) / (new_xwidth * new_height)

wszz = w_zfob(w13)/ (new_xwidth * new_ywidth)

wexx = 2.0 * xdif / (width + new_xwidth)

weyy = 2.0 * ydif / (width + new_ywidth)

wezz = 2.0 * zdif / (0.5*height + new_height)

wevol = wexx + weyy + wezz

end
; ---------------------------------------------



def make_gain ;

alpha = 0.05 ;

; count = 0

; avg_stiff = 0

; cp = contact_head ;

; loop while cp # null

; if c_ball2(cp) = w4

; count = count + 1

; avg_stiff = avg_stiff + c_kn(cp)

; end_if

; if c_ball2(cp) = w3

; count = count + 1

; avg_stiff = avg_stiff + c_kn(cp)

; end_if

; if c_ball2(cp) = w5

; count = count + 1

; avg_stiff = avg_stiff + c_kn(cp)

; end_if

; if c_ball2(cp) = w2

; count = count + 1

; avg_stiff = avg_stiff + c_kn(cp)

; end_if

; cp = c_next(cp)

; end_loop

;ncount = count / 4.0

;avg_stiff = avg_stiff / count

;gxy = alpha * (6.05 * 6.0) / (avg_stiff * ncount * tdel)

count = 0

avg_stiff = 0

cp = contact_head ; find avg. number of contacts on top/bottom walls

loop while cp # null

; if c_ball2(cp) = w1

; count = count + 1

; avg_stiff = avg_stiff + c_kn(cp)

; end_if

if c_ball2(cp) = w13

if wszz#0

www=c_ball2(cp)

count = count + 1

avg_stiff = avg_stiff + c_kn(cp)

command

;print count

;pause

end_command

end_if

end_if

cp = c_next(cp)

end_loop

ncount = count

if count # 0

avg_stiff = avg_stiff / count

gz = alpha * (6.05 * 6.0)/ (avg_stiff * ncount * tdel)

else

gz=0.1

end_if

command

print count gz

;pause

end_command

end
;-----------------------------------------------------------

def _servo

while_stepping

make_ss

ydisp=ydisp+w_y(w4)

; udx = gxy * (wsxx - 100000)

; w_xvel(w2) = udx

; w_xvel(w3) = -udx

; udy = gxy * (wsyy -  100000)

; w_yvel(w4) = udy

; w_yvel(w5) = -udy

; if z_servo = 1 ; switch stress servo on or off

  udz = gz * (wszz -  100000)

  command

  print udz gz wszz count

; pause

end_command

; w_zvel(w1) = udz

  w_zvel(w13) = udz

; end_if

end
;-----------------------------------------------------------

def iterate

loop while 1 # 0

n_cyc=n_cyc+1

make_gain

; if abs((wsxx -  100000)/ 100000) < 5 then

; if abs((wsyy - 100000)/ 100000) < 5 then

; if wszz> 100000

if abs((wszz -  100000)/ 100000) < 5 then

exit

end_if

; end_if

; end_if

; end_if

command

cyc 100

print udz gz wszz count

end_command

print udz gz wszz count

end_loop

end

;----------------------------------------------

def accel_platens

;------Accelerates the platens to achieve vel of _vfinal in _nsteps,

; using _nchunks

_nsteps=2000

_nchunks=200

_vfinal=0.02

_close=0

_niter=_nsteps/_nchunks

loop _chnk(1,_nchunks)

if _close=1 then

  _vel=_chnk*(_vfinal/_nchunks)

else

  _vel=-_chnk*(_vfinal/_nchunks)

end_if

  _mvel=-_vel

  command

  wall id 2 yv _vel

  wall id 4 yv _vel

  wall id 3 yv _vel
  
  wall id 5 yv _vel

  wall id 6 yv _vel

  wall id 7 yv _vel

  step _niter

end_command

end_loop

end

;----------------------------------------------------

set grav 0 0 -9.8

macro zero_vel 'ini xvel 0.0 yvel 0.0 zvel 0.0 xspin 0.0 yspin 0.0 zspin 0.0'

plot add axes brown

plot add wall white id on

;plot show


;zero_vel

;cyc 5

;zero_vel

;cyc 5

zero_vel

wall_addr


prop fric 0.7

wall_addr

ydisp1=w_y(w4)

ydisp2=w_y(w3)

xdif = w_x(w8) - w_x(w10)

ydif = w_y(w11)

zdif = w_z(w13)

wyforce2=w_yfob(w2)

wyforce8=w_yfob(w8)

wyforce10=w_yfob(w10)

wyforce11=w_yfob(w11)

wyforce9=w_yfob(w9)

wxforce2=w_xfob(w2)

wxforce8=w_xfob(w8)

wxforce10=w_xfob(w10)

wxforce11=w_xfob(w11)

wxforce9=w_xfob(w9)

shearforce=0.5*(wyforce2-wyforce11)

new_xwidth = 6.0 + xdif

new_ywidth = 6.0 + ydif

new_height = 0.5*6.0 + zdif

wsxx = 0.5*(w_xfob(w8) - w_xfob(w10)) / (new_ywidth * new_height)

wsyy = w_yfob(w11) / (new_xwidth * new_height)

wszz = w_zfob(w13)/ (new_xwidth * new_ywidth)

wexx = 2.0 * xdif / (width + new_xwidth)

weyy = 2.0 * ydif / (width + new_ywidth)

wezz = 2.0 * zdif / (0.5*height + new_height)

wevol = wexx + weyy + wezz

iterate ;

hist id=1 shearforce

hist id=2 wszz

hist id=3 ydisp2

hist id=4 wyforce2

hist id=5 wyforce11

hist id=6 wyforce8

hist id=7 wyforce9

hist id=8 wxforce2

hist id=9 wxforce8

hist id=10 energy body

hist id=11 energy boundary

hist id=12 energy fric

hist id=13 energy kinetic

hist id=14 energy strain

plot add hist -5 vs -3

accel_platens

save directshear_reg.sav