首先感谢中南王同学提供的滚刀图形,PFC外部导入wall,需要对网格有要求,需要都是三角形网格。所以画好了图形轮廓后需要使用犀牛软件调整一下网格才能够正确导入。
滚刀破岩应该算是这几年大家做的比较多的案例,其实这个案例的难点不算很多,第一个就是上面说的如何导入滚刀的wall,第二个就是滚刀滚动的额模拟。
案例的步骤分为:成样、预压、加胶结、加刀片、滚动。
前面三个和前面都一样,这里就不再赘述了。
看一下加刀片的代码:
首先指定了刀片的半径,这里用的也是14英寸刀片的半径,因为我们这里有删除墙的操作,所以需要将伺服函数移除,防止空指针报错。同时也不要忘了将墙的速度设为0。
restore jiajiaojie
[Rdao=0.432]
set fish callback -1.0 remove @sevro_walls
wall delete walls range id 2
wall delete walls range id 4
wall attribute velocity 0 0 0 range id 1
wall attribute velocity 0 0 0 range id 3
wall attribute velocity 0 0 0 range id 5
wall attribute velocity 0 0 0 range id 6
call addgeowall
@add_geo("daopan1.stl")
@zhihuanY_X(1)
[n_rad=Rdao/(z_max-z_min)]
@sacle(1,@n_rad,1)
@sacle(1,@n_rad,3)
@sacle(1,@n_rad,2)
@shengchengwall(1,[wlx*0.5],0,[wlz*0.5 Rdao*0.5],10)
save jiadaopian
这里引用addgeowall的是我之前写的一个dat文件,和geo_tools比较类似。首先引入daopian1的图形文件,因为画的时候没有注意坐标,所以导入的时候发现墙体方向不对,所以使用zhihuanY_X函数转了一下。之后就是设置一下放大系数合理的调整刀片的大小,并将其放置于试样的边角位置。代码如下:
def add_geo(ss)
command
geometry import @ss
endcommand
get_min_max(1)
moveToOrigin(1)
end
def get_min_max(id)
global x_min=1e100
global x_max=-1e100
global y_min=1e100
global y_max=-1e100
global z_min=1e100
global z_max=-1e100
local gs = geom.set.find(id)
loop foreach local gn geom.node.list(gs)
local pos = geom.node.pos(gn)
if x_min > comp.x(pos)
x_min = comp.x(pos)
endif
if y_min > comp.y(pos)
y_min = comp.y(pos)
endif
if z_min > comp.z(pos)
z_min = comp.z(pos)
endif
if x_max < comp.x(pos)
x_max = comp.x(pos)
endif
if y_max < comp.y(pos)
y_max = comp.y(pos)
endif
if z_max < comp.z(pos)
z_max = comp.z(pos)
endif
endloop
end
def moveToOrigin(id)
get_min_max(id)
local gs = geom.set.find(id)
loop foreach local gn geom.node.list(gs)
geom.node.pos.x(gn)=geom.node.pos.x(gn)-(x_max x_min)*0.5
geom.node.pos.y(gn)=geom.node.pos.y(gn)-(y_max y_min)*0.5
geom.node.pos.z(gn)=geom.node.pos.z(gn)-(z_max z_min)*0.5
endloop
get_min_max(id)
end
def sacle(id,n,option)
local gs = geom.set.find(id)
loop foreach local gn geom.node.list(gs)
if option=1 then
geom.node.pos.x(gn)=geom.node.pos.x(gn)*n
else if option=2 then
geom.node.pos.y(gn)=geom.node.pos.y(gn)*n
else
geom.node.pos.z(gn)=geom.node.pos.z(gn)*n
endif
endloop
get_min_max(1)
end
def move(id,x_mov,y_mov,z_mov)
local gs = geom.set.find(id)
loop foreach local gn geom.node.list(gs)
geom.node.pos.x(gn)=geom.node.pos.x(gn) x_mov
geom.node.pos.y(gn)=geom.node.pos.y(gn) y_mov
geom.node.pos.z(gn)=geom.node.pos.z(gn) z_mov
endloop
get_min_max(1)
end
def zhihuanY_Z(id)
local gs = geom.set.find(id)
loop foreach local gn geom.node.list(gs)
y_cord=geom.node.pos.y(gn)
z_cord=geom.node.pos.z(gn)
geom.node.pos.y(gn)=z_cord
geom.node.pos.z(gn)=y_cord
endloop
end
def zhihuanY_X(id)
local gs = geom.set.find(id)
loop foreach local gn geom.node.list(gs)
y_cord=geom.node.pos.y(gn)
x_cord=geom.node.pos.x(gn)
geom.node.pos.y(gn)=x_cord
geom.node.pos.x(gn)=y_cord
endloop
end
def shengchengwall(id,pos_x,pos_y,pos_z,wall_id)
move(1,pos_x,pos_y,pos_z)
local gs = geom.set.find(id)
geoname=geom.set.name(gs)
wall_name=string.build("ceqiang_%1",wall_id)
command
wall import geometry @geoname id @wall_id name @wall_name group ceqiang
endcommand
move(1,-1*pos_x,-1*pos_y,-1*pos_z)
end
处理完后的模型如图:
下面是滚动部分:
这里对应现实的刀片自转速度为6.28rad/s,公转线速度为0.04m/s,为了节省时间,扩大100倍。定义一下切割的距离,就能够算出计算的时间。这里需要注意的点是需要时刻更新刀片的滚动中心,不然刀片会飞掉。
restore jiadaopian
[n_kuoda=100]
[xuanzhuansulv=6.28*n_kuoda]
[qianjin=-0.04*n_kuoda]
[xiaya=qianjin*0.2]
[qianjinjuli=-length*0.8]
[final_time=qianjinjuli/qianjin]
ball attribute displacement multiply 0
set mech age 0
wall attribute velocity @qianjin 0 @xiaya yspin @xuanzhuansulv range id 10
cmat add 1 model linear method deformability emod 100e6 kratio 1.5 property fric 0.5 range contact type ball-facet
def daopian_wp
wpdaopian=wall.find(10)
end
@daopian_wp
def upadta_spin_center
x_center=wlx*0.5 qianjin*mech.age
z_center=wlz*0.5 Rdao*0.5 xiaya*mech.age
wall.rotation.center.x(wpdaopian)=x_center
wall.rotation.center.z(wpdaopian)=z_center
end
def cal_niuju
niuju=0
loop foreach ct wall.contactmap(wpdaopian)
forcrct=contact.force.global.x(ct)
liju=contact.pos.z(ct)-wall.pos.z(wpdaopian)
niuju =forcrct*liju
endloop
end
def jiance
forecDaoZ=wall.force.contact.z(wpdaopian)
forecDaoX=wall.force.contact.x(wpdaopian)
time_now=mech.age*1.0
weiyiX=qianjin*mech.age
weiyiZ=xiaya*mech.age
cal_niuju
end
set fish callback -1.0 @jiance
set fish callback -1.1 @upadta_spin_center
[baocunpinlv=final_time/20.0]
[time_record=mech.age-baocunpinlv-1]
[count=1]
def savefile
jiance
upadta_spin_center
if mech.age-time_record > baocunpinlv then
filename=string.build("jieguo%1",count)
command
save @filename
plot bitmap filename @filename
endcommand
time_record=mech.age
count =1
endif
end
set fish callback -1.0 @savefile
history delete
history id 1 @forecDaoZ
history id 2 @time_now
history id 3 @weiyiX
history id 4 @weiyiZ
history id 5 @forecDaoX
history id 6 @niuju
history nstep 50
call fracture.p3fis
@track_init
solve time @final_time
save resultTBD
我们看一下最终的破坏图:
当然做成动图也很好:
俯视图也能看出掘进过程:
当然力链图也能看出刀片的受力,受力主要在前半部分,还是比较合理的。
裂纹的产生也是我们需要分析的东西,这里将颗粒透明度调高:
下面看一下刀片的Z向受力图:
也可以看一下根据接触力计算的扭矩图:
这两个数值都是和实际相符的。这里就不去多做研究了。这里比较麻烦的一点是,对于刀片大小的改变,我们可以进行处理,但是如果对其形状进行改变的话,就需要画不同的刀片图了,建议或许可以将刀片简化一下,这样就可以研究一下变量的影响了。