【PFC6.0】随机多边形区域划分及颗粒填充

lobby 2022年9月9日浏览：1948 收藏：6

0 引言

在做岩石的时候，粗糙一点就认为岩石是一个均匀的整体，对整体参数进行标定，达到一个单元层次的宏观特性对应。精细点做的话就需要考虑矿物，这样我们会认为矿物内部是一个均匀的整体，这样的话就需要对岩石内部的区域进行划分，规定每种矿物成分的区域，然后用颗粒进行填充模拟。

本文主要目的是将区域进行随机多边形划分，并且往其中填充颗粒。使用到的技术主要是rblock中的merge命令。

1 区域划分

我们主要使用rblock的建模特性作为过渡，得到随机的多边形区域。

model new

model domain extent -10 10geometry set "box"geometry generate box -5 5

def par    origin_rad_min=0.2    origin_rad_max=origin_rad_min*1.5    scale=3    block_radius_max=origin_rad_max*scale    block_radius_min=origin_rad_min*scaleend@parrblock construct from-geometry "box"  minimum-edge  @origin_rad_min ...        maximum-edge @origin_rad_max group "origin" slot "1"

第一步是利用rblock construct命令可以对区域进行三角划分，如果可以直接对区域进行多边形划分肯定是更好的。但就我目前对rblock的理解来说，还做不到这种程度，目前只能做到对区域进行三角划分，这里可以指定三角形的边长范围。

我们生成了一堆三角形的rblock，如下图：

【PFC6.0】随机多边形区域划分及颗粒填充的图1

第二步就是用所得到的三角形进行随机的合并，这里我们用的是长方形区域进行合并。

def inbox(pos,pointer)    groupName=rblock.group(pointer)    if groupName#"1=origin" then        inbox=false        exit    endif    vl=vector(x_pos_random-sousuo_rad_x,y_pos_random-sousuo_rad_y)    vh=vector(x_pos_random+sousuo_rad_x,y_pos_random+sousuo_rad_y)    rbpos=rblock.pos(pointer)    if comp.x(rbpos)>comp.x(vl) & comp.x(rbpos)<comp.x(vh) &comp.y(rbpos)>comp.y(vl) & comp.y(rbpos)<comp.y(vh) then                inbox=true        exit    endif    inbox=falseend
def conbine_rblock    stop_flag=true    name_count=1    loop while stop_flag        x_pos_random=math.random.uniform*10-5        y_pos_random=math.random.uniform*10-5        sousuo_rad_x=math.random.uniform*(block_radius_max-block_radius_min)+block_radius_min        sousuo_rad_y=math.random.uniform*(block_radius_max-block_radius_min)+block_radius_min        ishaveRb=false        loop foreach rb rblock.list            vl=vector(x_pos_random-sousuo_rad_x,y_pos_random-sousuo_rad_y)            vh=vector(x_pos_random+sousuo_rad_x,y_pos_random+sousuo_rad_y)            if inbox(vl,rb) then                ishaveRb=true                exit loop            endif        endloop                group_name=string.build("zu_%1",name_count)        name_count+=1        if ishaveRb=true then                    command                                rblock merge group "new" slot "1" group @group_name slot "10"  ...                        range fish @inbox                model clean            endcommand        endif        orgin_count=0        loop foreach rb rblock.groupmap("origin",1)            orgin_count+=1        endloop        if orgin_count==0 then            stop_flag=false        endif    endloopend@conbine_rblock

我们每次随机获得一个长宽的长方形，将长方形中的“origin”组的rblock合并成一个新的rblock。

这里有个逻辑需要注意是，merge的概念不一定需要rblock紧贴在一起，就算rblock分开一段距离，他们之间的范围也会合并到新的rblock中，所以合并得到的rblock实际上是有重叠的。

我们看一下新的rblock。

【PFC6.0】随机多边形区域划分及颗粒填充的图2

可以发现的是，大部分的rblock都被合并了，但是有少部分rblock是被“剩余合并”的，所谓“剩余合并”的概念就是某个rblock或者某几个会被单独剩余下来，最终形成一个rblock，而这个rblock是不满足长宽要求的。

所以我们还需要进行第三步，也就是把这些小的rblock合并到它周围的大rblock中。

def findMinRb    VLimit=origin_rad_min*origin_rad_min*10    loop foreach rb rblock.groupmap("new",1)        if rblock.vol(rb)<VLimit*4 then            rblock.group(rb,1)="minRB"        endif    endloopend@findMinRb

def findNearestRb(pos)    minPos=1e100    minId=0    loop foreach rbfind rblock.groupmap("new",1)        vecLenghth=rblock.pos(rbfind)-pos        if math.mag(vecLenghth)<minPos then            minPos=math.mag(vecLenghth)            minId=rblock.id(rbfind)        endif    endloop    findNearestRb=minIdend

def solveMinRb        loop foreach rb  rblock.groupmap("minRB",1)        rbNearest=rblock.find(findNearestRb(rblock.pos(rb)))        rblock.group(rb,20)="dai"        rblock.group(rbNearest,20)="dai"        group_name=string.build("zu_%1",name_count)        name_count+=1        command             rblock merge group "new" slot "1" group @group_name ...                slot "10" range group "20=dai"        endcommand        loop foreach rbdai rblock.groupmap("dai",20)            rblock.group(rbdai,20)="none"        endloop    endloopend@solveMinRb

通过以上的命令，我们可以找到小的颗粒，当然我们也可以调整“小的颗粒”的范围，从而调整模型中区域的精细度。

运行完的效果如图：

【PFC6.0】随机多边形区域划分及颗粒填充的图3

这个就是我们区域的范围。

2 颗粒填充

由于没有指定随机数，导致本文出现的图和我既有算例的分布有区别，这里给出我算完的一个sav。如下图：

【PFC6.0】随机多边形区域划分及颗粒填充的图4

可以看到区域大小基本一致，只是形状不同，这里也体现随机性的概念。

model restore "rblock_slip"rblock export to-geometry "allgeo" cmat default type ball-facet model linear method deform emod 100e6 kratio 1.5 property fric 0.5cmat default type ball-ball model linear method deform emod 100e6 kratio 1.5  property fric 0.5def GetGeo    geoname_count=1    loop foreach rb rblock.list        groupName=string.build("geo_%1",geoname_count)        idRblock=rblock.id(rb)        command            rblock export to-geometry @groupName range id @idRblock        endcommand        geoname_count+=1    endloopend@GetGeorblock delete

颗粒填充前我们需要把rblock的边界都导出为geometry，方便后面根据区域进行颗粒生成。

可以看一下区域图：

【PFC6.0】随机多边形区域划分及颗粒填充的图5

由于rblock有重叠，导致geometry也是有重叠的，所以颗粒的生成规则我们需要进行优化。

首先定一个考虑区域重叠的规则。

区域交集优先被id小的geometry使用。

也就是说id为3区域需要减去其与id 1和id 2区域的交集，作为对应颗粒的生成范围。

榆树我们可以得到某个区域的颗粒投放逻辑，如下：

def ballDeleteGeo(geo_count_single)    if geo_count_single=1 then        exit    endif    thisGroup=string.build("geo_%1",geo_count_single)    loop n(1,geo_count_single-1)         groupName=string.build("geo_%1",n)        command            ball delete range geometry-space @groupName count 1 group @thisGroup        endcommand    endloopenddef GenerateBallIn(geo_count_single)    groupName=string.build("geo_%1",geo_count_single)    command        ball distribute porosity 0.1 radius 0.03 0.06 group @groupName range geometry-space @groupName count 1         ball attribute density 2.7e3 damp 0.3    endcommand    ballDeleteGeo(geo_count_single)end

第二个逻辑是颗粒的平衡边界，每个区域的颗粒投放后都需要平衡，其边界应当考虑id小于当前geo id的区域，所以需要控制墙体的生成和删除。

wall import  from-geometry "box" id 100 def genWalls(geo_count_single)    loop n(1,geo_count_single)         groupName=string.build("geo_%1",n)        command            wall import  from-geometry @groupName id [n+10000]         endcommand    endloopenddef deleWalls(geo_count_single)    loop n(1,geo_count_single)        command            wall delete walls range id [n+10000]         endcommand    endloopend

def genAllBall    loop gen_cur(1,geoname_count-1)                genWalls(gen_cur)        command                       model calm            ball fix velocity spin                   endcommand        GenerateBallIn(gen_cur)        groupName=string.build("geo_%1",gen_cur)        command                       model cycle 2000 calm 20            ball delete range geometry-space @groupName count 1 not group @groupName        endcommand        ballDeleteGeo(gen_cur)        command             model solve ratio-average 1e-5 or cycles 10000            model save @groupName        endcommand        deleWalls(gen_cur)    endloopend@genAllBallball free velocity spinmodel cycle 1model solvemodel save "sample"