PFC三维Cluster碎石三轴模拟实操（附赠程序代码）

详细信息

作者优秀: 优秀教师/意见领袖/博士学历/特邀专家/独家讲师
平台推荐: 内容稀缺

1年前浏览4565

导读：离散元作为一款基于颗粒力学的软件，得益于胶结模型的开发应用，使其在模拟胶结材料上展现出了令人惊讶的效果。对于散体材料，使用ball模拟可以很好的模拟颗粒在变形过程中的滑动、跃迁。但是美中不足的是，ball作为刚性体，不能体现出颗粒的破碎效果。

目前主流的是有两种方法，一种是监测ball的应力状态，当其超过强度阈值时，用若干小颗粒来代替大颗粒完成破碎效果。这种方法比较好的是计算效率高，且计算简单，毕竟都是线性模型。缺点也很明显，大颗粒碎成多少块小颗粒、每个小颗粒多大，这个事情讲不清楚的话肯定达不到一个比较好的模拟效果。

另一种就是用cluster来模拟了，所谓cluster，中文名为柔性簇。其作用方式是若干细小颗粒聚合成指定形状，且胶结在一起。这种方法模拟效果较好，且能够比较好的反映碎石破碎后的状态。

本文基于之前研究的成果 柔性簇Cluster模拟基本框架介绍 ，开发一个适用于三维碎石的cluster框架，并在三轴中进行实现。

一、框架介绍

本文和之前讲述的二维一样，首先进行模板生成，只是现在不是做cluster模板，而是做clump模板。完成一个比较好的clump模板后，在成样时候先完成clump的成样，再进行cluster颗粒的替换。

二、clump模板生成

clump中ball替换pebble的概念早已经有了，其实缺少的就是一个比较好的clump模板。clump中pebble不能重叠量过大，不然在替换成ball后ball也是有大的重叠量，而这种计算是不准确的。

所以这里我们必须想出一个生成pebble重叠量不大的clump模板。

最后采用的模拟思路是：1）将形状导入为墙体，并且在墙体中生成指定孔隙率的颗粒进行平衡；2）读取当前所有的ball作为clump模板；3）模板导出为本地的p3clp文件格式，后续用到成样工况中。

代码如下：

model newmodel domain extent -2 2def RdPar    rdmin=0.08    rdmax=0.12       poro=0.3end@RdParcmat default type ball-facet model linear property kn 1e8 cmat default type ball-ball model linear property kn 1e7 ks 1e7 fric 0.2 program call "Geo_Tools.dat"program call "OutClumpTemplate.dat"def CreateTemplate    loop n(1,4)        stlName="Rock"+string(n)        fileName="shape\\"+stlName+".stl"        command            geometry import @fileName                    endcommand
        get_min_max(stlName)        x_len=x_max-x_min        move_geo(stlName,vector(-(x_max+x_min)*0.5,-(y_max+y_min)*0.5,-(z_max+z_min)*0.5))        scale_geo(stlName,1.0/x_len)        command            wall import  from-geometry @stlName            ball distribute radius @rdmin @rdmax porosity @poro range geometry-space @stlName inside            ball attribute density 2e3 damp 0.7            model cycle  2000 calm 50            model cycle  10000        endcommand        OutTemplate(stlName)        command            model save @stlName            ball delete            wall delete            clump template delete        endcommand    endloopend@CreateTemplate

其中先假设我们的形状文件存在于本地的shape文件夹中，且以Rock来命名，这里可以自行修改。

文中用到的Geo_tools文件内容如下：

主要的目标是实现导入的图形缩放到x向长度为1，会在导入后先move到原点，再进行scale.

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    endloopend

def scale_geo(geo_name,scale_factor)    geo_zidan=geo.set.find(geo_name)    loop foreach nd geo.node.list(geo_zidan)        geo.node.pos.x(nd)=geo.node.pos.x(nd)*scale_factor        geo.node.pos.y(nd)=geo.node.pos.y(nd)*scale_factor        geo.node.pos.z(nd)=geo.node.pos.z(nd)*scale_factor    endloopend

def move_geo(geo_name,dist)    geo_zidan=geo.set.find(geo_name)    loop foreach nd geo.node.list(geo_zidan)        geo.node.pos(nd)=geo.node.pos(nd)+dist    endloopend

文中用到的第二个文件OutClumpTemplate代码如下:

这里的目的是读取当前所有的颗粒，在当前目录的Template文件夹中生成clump模板的p3clp后缀的文件。

def OutTemplate(templateName)    ball_num=ball.num    command        clump template create geometry @templateName  bubblepack ratio 0.6 distance 50 ...                surfcalculate    endcommand
    cp_template= clump.template.find(templateName)    loop foreach pb  clump.template.pebblelist(cp_template)        clump.template.deletepebble(cp_template,pb)    endloop    loop foreach bp ball.list        clump.template.addpebble(cp_template,ball.radius(bp),ball.pos(bp))    endloop       clump_file_name="Template\\"+templateName+".p3clp"    command        clump template export @templateName filename @clump_file_name    endcommandend

到这一步，我们的四个模板就生成了。效果如下：

Rock1:

Rock2:

Rock3:

Rock4:

以上就是达到我们使用的clump模板要求，即pebble不可有大的重叠量。

p3clp文件存储路径如下

然后下面进行三轴的时候，把这个文件夹移动到三轴项目目录中

三、碎石三轴----成样

生成指定粒径、空隙率、尺寸的式样，这部分略过了。

model newdef par    width=2    height=width*2
    rdmin=0.15    rdmax=0.2
    poro=0.16end@par

domain extent [-width*2.0] [width*2.0] [-width*2.0] [width*2.0] ...        [-height*2.0] [height*2.0]model random 10001

wall generate box [-width*0.5] [width*0.5] [-width*0.5] [width*0.5] ...            [-height*0.5] [height*0.5] expand 1.5ball distribute porosity @poro radius [rdmin] [rdmax] box ...    [-width*0.5] [width*0.5] [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] ball attribute density 2e3 damp 0.7contact cmat default model linear method deform emod 100e6 kratio 1.5model cycle 2000 calm 50ball property "fric" 0.5model solvemodel save "sample"

四、碎石三轴----clump替换

先进行第一种替换：ball替换为clump，这里其实也是用clump作为桥梁来确定cluster中颗粒的位置。

这部分代码和之前三点弯曲基本类似。

有一点是开头首先进行clump template import，导入本地的模板文件。

还有就是模板选择，这里有四个模板，目前完成的概念是四个模板随机均匀分布，如果有规则的话，可以修改0.25 0.5 0.75这三个界限值。

注意避免墙体和clump自锁，所以墙体先放大，再压缩。

model restore "sample"

def importTemplate    loop n(1,4)        stlName="Rock"+string(n)        fileName="Template\\"+stlName+".p3clp"        command            clump template import @fileName name @stlName        endcommand    endloopend@importTemplate

def GetTemplateName    rd=math.random.uniform    if rd<0.25 then        GetTemplateName="Rock1"    else if rd<0.5 then         GetTemplateName="Rock2"    else if rd<0.75 then         GetTemplateName="Rock3"    else         GetTemplateName="Rock4"    endifend

[clump_count=1]def tihuan    loop foreach bp ball.list        x_pos   = ball.pos(bp,1)        y_pos   = ball.pos(bp,2)        z_pos   = ball.pos(bp,3)        bvol = (4/3.0)*math.pi*ball.radius(bp)^3           template_name=GetTemplateName        ball.delete(bp)  
        angle= math.random.uniform*180         axis_x=math.random.uniform-0.5        axis_y=math.random.uniform-0.5        axis_z=math.random.uniform-0.5        axis=vector(axis_x,axis_y,axis_z)        
        command                            clump replicate id @clump_count name @template_name  ...                            pos-x @x_pos pos-y @y_pos pos-z @z_pos  ...                           volume @bvol  angle @angle axis @axis            clump group  @template_name range id @clump_count        endcommand          clump_count+=1    endloopend@tihuan

clump attribute density 2e3 damp 0.7

cmat default type pebble-facet model linear method deform emod 1000e6 kratio 1.5clump attribute density 2.3e3 damp 0.7





[yasuo_time=5]

wall deletewall generate box [-width] [width] [-width] [width] ...    [-height] [height]  wall attribute velocity-z [(height*0.5)/yasuo_time] range id 1wall attribute velocity-z [-(height*0.5)/yasuo_time] range id 2wall attribute velocity-x [(width*0.5)/yasuo_time] range id 3wall attribute velocity-x [-(width*0.5)/yasuo_time] range id 4 wall attribute velocity-y [(width*0.5)/yasuo_time] range id 5wall attribute velocity-y [-width*0.5/yasuo_time] range id 6solve time [yasuo_time*0.3] calm 10 solve time [yasuo_time*0.7]wall attribute velocity 0 0 0

model solve 

model save "tihuan_clump"

五、碎石三轴----cluster替换

这里是我们主要的工作量了，模拟的概念是：

1）找到每一个clump，找到这个clump中每一个pebble的位置和半径，在这个位置上生成同粒径的ball

2）把这些颗粒打个组叫“jiaojie”，然后即时性的给这个组的颗粒附上pb模型，并且加上接触。由于只针对这个组，且指定了match 2，所以“jiaojie”这个组和其余的组之间的接触走default生成linear接触。每次进行内应力的清零防止颗粒崩散。安全起见也fix一下。

3）这些颗粒之间的接触数值可在外面写函数获取，比如我这里的GetEmodByClusterName,这里传进来的参数是cluster的group，也就是模板的名称，这里认为同一个模板用的接触属性一样，可以自己去确定属性规则。

model restore "tihuan_clump"

contact cmat default type ball-ball model linear method deform emod 1e9 kratio 1.5 property fric 0.2 lin_mode 1contact cmat default type ball-facet model linear method deform emod 10e9 kratio 1.5 property fric 0.2 lin_mode 1def GetEmodByClusterName(clusterName)    test_string=clusterName    if clusterName=="Rock1" then        GetEmodByClusterName=1e9    else if clusterName=="Rock2" then        GetEmodByClusterName=2e9    else if clusterName=="Rock3" then        GetEmodByClusterName=3e9    else        GetEmodByClusterName=4e9    endifend

def GetPbCohByClusterName(clusterName)    if clusterName=="Rock1" then        GetPbCohByClusterName=1e7    else if clusterName=="Rock2" then        GetPbCohByClusterName=2e7    else if clusterName=="Rock3" then        GetPbCohByClusterName=3e7    else        GetPbCohByClusterName=4e7    endifend

def applyBond(groupName,clusterName,bondGap)    emod=GetEmodByClusterName(clusterName)    pb_coh=GetPbCohByClusterName(clusterName)    command        model cycle 1        model calm        ball attribute force-contact multiply 0.0         clump attribute force-contact multiply 0.0         contact property lin_force 0.0 0.0 0.0 range group @groupName        contact model linearpbond range group @groupName match 2        contact method deform emod [emod] kratio 1.5 pb_deform emod [emod] kratio 1.5 range group @groupName match 2        contact property pb_coh @pb_coh pb_ten [pb_coh/1.5] pb_fa 50 fric 0.1 lin_mode 1 range group @groupName match 2        model clean        contact method bond gap @bondGap range group @groupName match 2    endcommandend



def TihuanCp    loop foreach cp clump.list        templateName=clump.template.name(clump.template(cp))         min_rd=1e100        loop foreach pb clump.pebblelist(cp)            rd=clump.pebble.radius(pb)            min_rd=math.min(min_rd,rd)            pos=clump.pebble.pos(pb)            bp=ball.create(rd,pos)            ball.group(bp)="jiaojie"        endloop        command            ball attribute density 2e3 damp 0.7 range group "jiaojie"        endcommand        applyBond("jiaojie",templateName,min_rd)        clump.delete(cp)        command            ball fix vel spin range group "jiaojie"            ball group @templateName range group "jiaojie"        endcommand    endloopend@TihuanCp

model save "tihuan_cluster"

生成完后，显示接触模型名称进行查看：

可以明显看到只有碎石内部是有胶结的。

显示pb_ten:

可以看到不同的属性也是给进去的。

六、碎石三轴----加围压

没什么可讲的，简单的伺服概念

model restore "tihuan_cluster"ball free vel spin [txx=-2e6][tyy=-2e6][tzz=-2e6]

[sevro_factor=0.5][do_xSevro=true][do_zSevro=true]

[sevro_freq=100][timestepNow=global.step-1]def sevro_walls    compute_stress     if timestepNow<global.step then        get_g(sevro_factor)        timestepNow+=sevro_freq    endif    if do_xSevro=true then        Xvel=gx*(wxss-txx)        Yvel=gy*(wyss-tyy)        vel_arg=(Xvel+Yvel)*0.5        wall.vel.x(wpRight)=-vel_arg        wall.vel.x(wpLeft)=vel_arg
        wall.vel.y(wpFront)=vel_arg        wall.vel.y(wpBack)=-vel_arg    endif
    if do_zSevro=true then        Zvel=gz*(wzss-tzz)        wall.vel.z(wpUp)=-Zvel        wall.vel.z(wpDown)=Zvel    endif
end

def wp_ini    wpDown=wall.find(1)    wpUp=wall.find(2)    wpLeft=wall.find(3)    wpRight=wall.find(4)    wpFront=wall.find(5)    wpBack=wall.find(6)end@wp_ini

def computer_chiCun    wlx=wall.pos.x(wpRight)-wall.pos.x(wpLeft)    wlz=wall.pos.z(wpUp)-wall.pos.z(wpDown)    wly=-wall.pos.y(wpFront)+wall.pos.y(wpBack)end

def compute_stress    computer_chiCun    wxss=-(wall.force.contact.x(wpRight)-wall.force.contact.x(wpLeft))*0.5/(wly*wlz)    wzss=-(wall.force.contact.z(wpUp)-wall.force.contact.z(wpDown))*0.5/(wlx*wly)    wyss=(wall.force.contact.y(wpFront)-wall.force.contact.y(wpBack))*0.5/(wlx*wlz)end

def get_g(fac)    gx=0    gy=0    gz=0    gX=1e10    gY=1e10    gZ=1e10    loop foreach ct wall.contactmap(wpLeft)        gX+=contact.prop(ct,"kn")    endloop    loop foreach ct wall.contactmap(wpRight)        gX+=contact.prop(ct,"kn")    endloop    loop foreach ct wall.contactmap(wpUp)        gZ+=contact.prop(ct,"kn")    endloop    loop foreach ct wall.contactmap(wpDown)        gZ+=contact.prop(ct,"kn")    endloop    loop foreach ct wall.contactmap(wpFront)        gY+=contact.prop(ct,"kn")    endloop    loop foreach ct wall.contactmap(wpBack)        gY+=contact.prop(ct,"kn")    endloop
    gx=fac*wly*wlz/(gX*global.timestep)    gy=fac*wlx*wlz/(gY*global.timestep)    gz=fac*wlx*wly/(gZ*global.timestep)end@get_g(@sevro_factor)

@compute_stress

fish callback add @sevro_walls -1.0 history deletehistory id 1 @wxsshistory id 2 @wysshistory id 3 @wzss;model mechanical timestep fix 1e-2

model cycle 1000

model solvemodel save "weiya"

七、碎石三轴----加载

model restore "weiya"

ball attribute displacement multiply 0 

[IZ0=wall.pos.z(wpUp)-wall.pos.z(wpDown)]

[strainrate=0.1][do_zSevro=false]wall attribute vel-z [-IZ0*strainrate*0.5] range id 2wall attribute vel-z [IZ0*strainrate*0.5] range id 1

def jiance    whilestepping    wezz=((wall.pos.z(wpUp)-wall.pos.z(wpDown))-IZ0)/IZ0end



history deletehistory id 1 @wzsshistory id 2 @wezz[final_time=0.1/strainrate]

[baocunpinlv=final_time/20.0][time_record_sav=-100][count=0]def savefile    time=mech.time.total    if time-time_record_sav >= baocunpinlv then        filename=string.build("jieguo_%1",count)        command           model save @filename        endcommand        time_record_sav=time        count +=1    endifend

fish callback add @savefile -1.0

program  call "fracture.p3fis"@track_init

model solve time [final_time]

model save "result"

首先看一下应力应变曲线：

比较经典的脆性破坏曲线。

所有颗粒在压缩过程中的破碎情况。

Rock1到Rock4的参数是逐渐变强的，所以碎石的破碎情况也是不均匀的。利于plot中的range筛选工具，可以选择只看某一个分组的变化情况。

Rock1

Rock2

Rock3

Rock4

可以看出Rock1由于参数最弱，所以坏的最多。

我这里裂纹文件没有改好，导致裂纹位置没有实现更新，就不分享出来了，大家可以自行用example中的裂纹文件。

从颗粒位置和裂纹位置的重合度来看，Rock1也是破碎最多的。

大家还可以去监测四个分组的破碎率，即胶结破坏数/总胶结数。计算不复杂，当作各位的作业了。

比较难的点在于级配统计了，这个也是一个比较广泛的需求。压缩破坏前后的级配曲线是我们的研究重点。

我这里也给出了计算级配的方法：

逻辑不复杂，用fragment的体积计算等体积下的球形颗粒半径作为粒径。然后统计当前所有fragment粒径范围，再进行分割统计。

model restore "jieguo_2"[split_num=20]

def GetFragmentDia(fg)    vol=fragment.vol(fg)    GetFragmentDia= (vol*3.0/4.0/math.pi)^(1/3.0)end



def GetMinMaxDiam    allVol=0    min_d=1e100    max_d=-1e100    loop foreach fg fragment.map        fg_vol=GetFragmentDia(fg)        allVol+=fragment.vol(fg)        min_d=math.min(min_d,fg_vol)        max_d=math.max(max_d,fg_vol)    endloop   end@GetMinMaxDiam[split_d=(max_d-min_d)/float(split_num)]def CreateTable    tb=table.create("jipei")    loop n(1,split_num)        cur_d=min_d+n*split_d        table.x(tb,n)=cur_d        table.y(tb,n)=0    endloopend@CreateTable

def tongji    loop foreach fg fragment.map        fg_vol=GetFragmentDia(fg)        idx=math.ceiling((fg_vol-min_d)/split_d)        if idx==0 then            idx=1        endif        loop n(idx,split_num)            y_value=table.y(tb,n)            y_value+=(fragment.vol(fg)/allVol)            table.y(tb,n)=y_value        endloop    endloop    end@tongji