本人在做一个拉拔试验模拟的时候发现位移和拉拔力始终不对，希望高手指点，，附上命令流
n
impgrid 3g
;==================================
def xiugai_11
p_gp=gp_head
loop while p_gp # null
w_temp1=gp_xpos(p_gp)
w_temp2=gp_ypos(p_gp)
gp_xpos(p_gp)=gp_zpos(p_gp)
gp_ypos(p_gp)=w_temp1
gp_zpos(p_gp)=w_temp2
p_gp = gp_next(p_gp)
endloop xiugai_1
end
xiugai_11
;==================================
gen separate 1
ini x add 10 range gro 1
int 1 face range cyl end1 (0 0 0) end2 (5 0 0) radius 0.00763
pl inter
pl bl gr
ini x add 10 range gro 1
ini x add -20 range gro 1
int 1 wrap 1 2
inter 1 maxedge 0.3
int 1 prop ks 1e10 kn 1e10 coh 1080 fric 10
;===================================
m m
;ini dens 7800 range gro 1
;ini dens 1290 range gro 2
;ini dens 1490 range gro 3
;ini dens 2060 ran gro 4
;ini dens 1930 ran gro 5
;set grav 0 0 -9.81
pro bulk 1.515e8 shear 7.8e7 tens 1.82e6 range gro 1
pro bulk 2.65e6 shear 1.747e6 tens 9e2 range gro 2
pro bulk 1.67e6 shear 7.7e5 tens 1.3e8 range gro 3
pro bulk 6.25e6 shear 5.08e6 tens 4e3 range gro 4
pro bulk 66.7e3 shear 30.7e3 tens 1.5e2 coh 27 fri 25 range gro 5
;============================================
fix x range x -0.1 0.1
;fix y z range gro 1
;===============================
def initval
pull_vel=-1.0e-6
my_sum=0.0
end
;===============================
initval
ini xvel=pull_vel range gro 1 x -0.1 0.1
;====================================================================
def force_len
whilestepping
gpp=gp_head
loop while gpp # null
w_tol=(sqrt(gp_ypos(gpp)*gp_ypos(gpp)+gp_zpos(gpp)*gp_zpos(gpp))-0.00763);
if w_tol<=0 then
_sum=_sum+gp_xfunbal(gpp)
end_if
gpp=gp_next(gpp)
end_loop
force_len=_sum
disp_end=step*(-1.0)*pull_vel
end
;====================================================================
history nstep=10
hist unbal
hist force_len
hist gp xdis 0 0 0
pl hist 2 vs 3
step 100000
