首先感谢中南王同学提供的滚刀图形,PFC外部导入wall,需要对网格有要求,需要都是三角形网格。所以画好了图形轮廓后需要使用犀牛软件调整一下网格才能够正确导入。
滚刀破岩应该算是这几年大家做的比较多的案例,其实这个案例的难点不算很多,第一个就是上面说的如何导入滚刀的wall,第二个就是滚刀滚动的额模拟。
案例的步骤分为:成样、预压、加胶结、加刀片、滚动。
前面三个和前面都一样,这里就不再赘述了。
看一下加刀片的代码:
首先指定了刀片的半径,这里用的也是14英寸刀片的半径,因为我们这里有删除墙的操作,所以需要将伺服函数移除,防止空指针报错。同时也不要忘了将墙的速度设为0。
restore jiajiaojie[Rdao=0.432]set fish callback -1.0 remove @sevro_wallswall delete walls range id 2wall delete walls range id 4wall attribute velocity 0 0 0 range id 1wall attribute velocity 0 0 0 range id 3wall attribute velocity 0 0 0 range id 5wall attribute velocity 0 0 0 range id 6call 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)commandgeometry import @ssendcommandget_min_max(1)moveToOrigin(1)enddef get_min_max(id)global x_min=1e100global x_max=-1e100global y_min=1e100global y_max=-1e100global z_min=1e100global z_max=-1e100local 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)endifif y_min > comp.y(pos)y_min = comp.y(pos)endifif z_min > comp.z(pos)z_min = comp.z(pos)endifif x_max < comp.x(pos)x_max = comp.x(pos)endifif y_max < comp.y(pos)y_max = comp.y(pos)endifif z_max < comp.z(pos)z_max = comp.z(pos)endifendloopenddef 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.5geom.node.pos.y(gn)=geom.node.pos.y(gn)-(y_max y_min)*0.5geom.node.pos.z(gn)=geom.node.pos.z(gn)-(z_max z_min)*0.5endloopget_min_max(id)enddef sacle(id,n,option)local gs = geom.set.find(id)loop foreach local gn geom.node.list(gs)if option=1 thengeom.node.pos.x(gn)=geom.node.pos.x(gn)*nelse if option=2 thengeom.node.pos.y(gn)=geom.node.pos.y(gn)*nelsegeom.node.pos.z(gn)=geom.node.pos.z(gn)*nendifendloopget_min_max(1)enddef 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_movgeom.node.pos.y(gn)=geom.node.pos.y(gn) y_movgeom.node.pos.z(gn)=geom.node.pos.z(gn) z_movendloopget_min_max(1)enddef 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_cordgeom.node.pos.z(gn)=y_cordendloopenddef 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_cordgeom.node.pos.x(gn)=y_cordendloopenddef 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)commandwall import geometry @geoname id @wall_id name @wall_name group ceqiangendcommandmove(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 0set mech age 0wall attribute velocity @qianjin 0 @xiaya yspin @xuanzhuansulv range id 10cmat add 1 model linear method deformability emod 100e6 kratio 1.5 property fric 0.5 range contact type ball-facetdef daopian_wpwpdaopian=wall.find(10)end@daopian_wpdef upadta_spin_centerx_center=wlx*0.5 qianjin*mech.agez_center=wlz*0.5 Rdao*0.5 xiaya*mech.agewall.rotation.center.x(wpdaopian)=x_centerwall.rotation.center.z(wpdaopian)=z_centerenddef cal_niujuniuju=0loop foreach ct wall.contactmap(wpdaopian)forcrct=contact.force.global.x(ct)liju=contact.pos.z(ct)-wall.pos.z(wpdaopian)niuju =forcrct*lijuendloopenddef jianceforecDaoZ=wall.force.contact.z(wpdaopian)forecDaoX=wall.force.contact.x(wpdaopian)time_now=mech.age*1.0weiyiX=qianjin*mech.ageweiyiZ=xiaya*mech.agecal_niujuendset fish callback -1.0 @jianceset fish callback -1.1 @upadta_spin_center[baocunpinlv=final_time/20.0][time_record=mech.age-baocunpinlv-1][count=1]def savefilejianceupadta_spin_centerif mech.age-time_record > baocunpinlv thenfilename=string.build("jieguo%1",count)commandsave @filenameplot bitmap filename @filenameendcommandtime_record=mech.agecount =1endifendset fish callback -1.0 @savefilehistory deletehistory id 1 @forecDaoZhistory id 2 @time_nowhistory id 3 @weiyiXhistory id 4 @weiyiZhistory id 5 @forecDaoXhistory id 6 @niujuhistory nstep 50call fracture.p3fis@track_initsolve time @final_timesave resultTBD
我们看一下最终的破坏图:

当然做成动图也很好:
.gif?imageView2/1/w/748/h/580)
俯视图也能看出掘进过程:
.gif?imageView2/1/w/759/h/522)
当然力链图也能看出刀片的受力,受力主要在前半部分,还是比较合理的。
.gif?imageView2/1/w/758/h/577)
裂纹的产生也是我们需要分析的东西,这里将颗粒透明度调高:
.gif?imageView2/1/w/732/h/592)
下面看一下刀片的Z向受力图:

也可以看一下根据接触力计算的扭矩图:

这两个数值都是和实际相符的。这里就不去多做研究了。这里比较麻烦的一点是,对于刀片大小的改变,我们可以进行处理,但是如果对其形状进行改变的话,就需要画不同的刀片图了,建议或许可以将刀片简化一下,这样就可以研究一下变量的影响了。