首页/文章/ 详情

PFC三维滚刀模拟

lobby

2年前浏览4600

首先感谢中南王同学提供的滚刀图形，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 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)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_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)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_cord        geom.node.pos.z(gn)=y_cord    endloopend

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    endloopend

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 0wall 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    endloopend
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_niujuend
set 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 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    endifend
set fish callback -1.0 @savefile
history deletehistory id 1 @forecDaoZhistory id 2 @time_nowhistory id 3 @weiyiXhistory id 4 @weiyiZhistory id 5 @forecDaoXhistory id 6 @niuju
history nstep 50call fracture.p3fis@track_init
solve time @final_time
save resultTBD