请教各位帮忙看看，我用的是颗粒球的加载，不是伺服，因为自己写的一个墙体伺服加载错了运行不了，就改用球的了（希望好心人有墙体加载的直剪命令发我下哈），出现了下面两个错误。
1.直剪的结果是上盒的球体往上走了。
2.上盒球体加的力该怎么算成应力。prop yforce=-1e6 range x 0,0.07 y 0.059,0.06
（添加不来附件啊，不好意思了，复制到下面）
模型建立
;
;fname: biax_1.DAT Creation and packing of particles - biaxial sample
new
set random ;reset random-number generator
set disk on ;treat balls as disks of unit thickness
; ;-----------------------------------------------------------------------------------------------------
;creat walls with overhang of extend
def make_walls
extend = 1.0
_x0 = width
_y0 = height/2
_x1 = width
_y1 = height* (1.0 + extend)
command
wall id = 1 kn = w_stiff nodes (_x0,_y0) (_x1,_y1)
end_command
_x0 = -width/2
_y0 = 0.0
_x1 = width* (1.0 + extend/2)
_y1 = 0
command
wall id = 2 kn = w_stiff nodes (_x0,_y0) (_x1,_y1)
end_command
_x0 = 0.0
_y0 = height * (1.0 + extend)
_x1 = 0.0
_y1 = height/2
command
wall id = 3 kn = w_stiff nodes (_x0,_y0) (_x1,_y1)
end_command
_x0 = width* (1.0 + extend/2)
_y0 = height
_x1 = -width/2
_y1 = height
command
wall id = 4 kn = w_stiff nodes (_x0,_y0) (_x1,_y1)
end_command
_x0 = 0.0
_y0 = height/2
_x1 = 0.0
_y1 = -height
command
wall id = 5 kn = w_stiff nodes (_x0,_y0) (_x1,_y1)
end_command
_x0 = width
_y0 = -height
_x1 = width
_y1 = height/2
command
wall id = 6 kn = w_stiff nodes (_x0,_y0) (_x1,_y1)
end_command
_x0 = -width
_y0 = height/2
_x1 = 0
_y1 = height/2
command
wall id = 7 kn = w_stiff nodes (_x0,_y0) (_x1,_y1)
end_command
_x0 = width
_y0 = height/2
_x1 = width* (1.0 + extend)
_y1 = height/2
command
wall id = 8 kn = w_stiff nodes (_x0,_y0) (_x1,_y1)
end_command
end
; ;-----------------------------------------------------------------------------------------------------
; assemble sample
def assemble
s_stiff = 0.0 ;shear stiffness
n_stiff = 5e8 ;normal contact stiffness
w_stiff = 5e9 ;
tot_vol = height * width * 1.0
rbar=0.5 * (rlo + rhi) ;求颗粒的平均半径
num = int((1.0 - poros) * tot_vol / (pi*rbar^2)) ;求可生成颗粒的总数目
mult = 1.6 ;initial radius multiplication factor
rlo_0 = rlo / mult ;initial minimum radius
rhi_0 = rhi / mult ;initial maximum radius
make_walls
command
gen id=1,num rad=rlo_0,rhi_0 x=0,width y=0,height tries 10000000
prop dens=1000 ks=s_stiff kn=n_stiff
end_command
ii = out(string(num)+' particles were created')
; get actual porosity
sum = 0.0
bp = ball_head
loop while bp # null ;通过遍历所有颗粒计算颗粒总体积
sum = sum + pi * b_rad(bp)^2
bp = b_next(bp)
end_loop
pmeas = 1.0 - sum / tot_vol ;actual porosity
mult = sqrt((1.0 - poros) / (1.0 - pmeas))
command
ini rad mul mult
cycle 1000
prop ks=5e8 fric 0.25
cycle 250
end_command
end
set dt dscale
print info
; ;-----------------------------------------------------------------------------------------------------
; change lateral wall stiffnesses
def cws
command
wall id 2 kn=w_stiff
wall id 4 kn=w_stiff
end_command
end
; ;-----------------------------------------------------------------------------------------------------
macro zero 'ini xvel 0 yvel 0 spin 0' ;通过采用宏macro将命令'ini xvel 0 yvel 0 spin 0'简化为zero
SET height=0.06 width=0.07 rlo=0.0003 rhi=0.0004 poros=0.1
assemble
SET w_stiff= 5e8 ; make lateral wall stiffness=1/10 of ball stiffness
cws
cyc 500
zero
; ;-----------------------------------------------------------------------------------------------------
def bound
bp=ball_head
loop while bp#null ;scan all balls
section
cp=b_clist(bp)
loop while cp # null ;scan ball's contacts
if c_nforce(cp) # 0.0 then
b2=c_ball2(cp)
if pointer_type(b2)=101 then ;b2 is a wall
b_xfix(bp)=1 ;fix original ball in x,y
b_yfix(bp)=1
b_color(bp)=1
exit section ;all done for this ball
end_if
end_if
if c_ball1(cp)=bp
cp=c_b1clist(cp)
else
cp=c_b2clist(cp)
end_if
end_loop
end_section
bp=b_next(bp)
end_loop
end
; ;-----------------------------------------------------------------------------------------------------
bound
ini xvel=0 yvel=0 spin=0
del wall 1
del wall 2
del wall 3
del wall 4
cycle 1
; ;-----------------------------------------------------------------------------------------------------
def replace ;Replace unbalanced forces wiht applied forces
bp=ball_head
loop while bp # null
if b_xfix(bp)=1
b_xfap(bp)=-b_xfob(bp)
b_yfap(bp)=-b_yfob(bp)
end_if
bp=b_next(bp)
end_loop
end
;
def qingling
bp=ball_head
loop while bp # null
b_xvel(bp)=0
b_yvel(bp)=0
b_xdisp(bp)=0
b_ydisp(bp)=0
bp=b_next(bp)
end_loop
end;-----------------------------------------------------------------------------------------------------
JSET id=1 dip=0 origin=(0.035,0.03)
replace
qingling
save MX.SAV
pause
; ;-----------------------------------------------------------------------------------------------------
; ;给定力学参数
prop fric=0.5
prop pb_kn=1e8 pb_ks=1e8 pb_rad=0.5
prop pb_nstren=5e7 pb_sstren=5e7
property n_bond=0 s_bond=0 range jset 1
property fric 0.2 range jset 1
plot creat bound_condition
plot add ball blue red
;plot add vel black
;plot add cf white
plot show
cycle 1000
; ;-----------------------------------------------------------------------------------------------------
;-----------------------------------------------------------------------------------------------------
; ;施加边界条件
free x y range x 0.001,0.069 ;leave ends fixed
ini xvel=1e-6 range x 0,0.001 y 0.03,0.06
fix x range x 0,0.001 y 0.03,0.06
ini xvel=1e-6 range x 0,0.07 y 0.059,0.06
fix x range x 0,0.07 y 0.059,0.06
ini xvel=1e-6 range x 0.069,0.07 y 0.03,0.06
fix x range x 0.069,0.07 y 0.03,0.06
prop yforce=-1e6 range x 0,0.07 y 0.059,0.06
fix y range x 0,0.07 y 0,0.001
; ;-----------------------------------------------------------------------------------------------------
cycle 1000
