简述PFC成样方法

    这部分内容包括：

    粒径放大法、压缩法、分层压缩、distribute成样法、分层成样法、Brick方法

    PFC生成颗粒的方法其实经历过了很长时间的发展，这里从发展的逻辑顺序来讲一下PFC中的成样方法。虽然是简述，但是这部分内容比较多，包含了我阅读的文献内容和自己经验和理解，建议或许可以打印下来学习理解一下。

    密样因为接触力的调整，最后肯定会比较均匀，所以松样才是衡量成样好坏的标准。

    我们这里生成的试样为：

    长宽为1.5m，两个粒径的颗粒，一个是0.006半径的，一个是0.009半径的，两个体积百分比一样，孔隙率为0.25。

一、generate方法

    generate方法为在指定区域生成指定数量的无重叠的颗粒（这句话理解一下）。可以看做是往一个棋盘上摆棋子这种感觉，颗粒是一个个生成的。逻辑方法为：1）生成随机数，确定区域内的坐标；2）判断坐标位置生成颗粒的话，会不会和别的颗粒重叠，如果重叠，回到第一步。

    这里就会暴露出generate方法的一个问题，就是可能在指定的循环次数内，找不到符合条件的坐标了，这样在区域内生成的颗粒就会远远少于指定的数目，下面用个demo说明一下这个问题。

    根据试样的条件，我们可以计算出两个粒径的颗粒数。逻辑为：试样颗粒体积 Vs=(1-poro)*V ，其中poro为孔隙率，V为试样体积。如果两个颗粒的质量百分比一样的话（如果密度一样，体积百分比也一样），两个粒径的颗粒体积为 Vs1=Vs2=Vs*0.5 。一个颗粒的体积为 pi*r*r ，这样颗粒数为 Vs1/(pi*r*r)。于是乎命令流为：

new def par    width=1.5    height= width       rdmax=0.009    rdmin=0.006        poro=0.25    rdmax_num=math.floor((1-poro)*0.5*width*height/(math.pi*rdmax*rdmax))    rdmin_num=math.floor((1-poro)*0.5*width*height/(math.pi*rdmin*rdmin))end@pardomain extent [-height*5] [height*5]set random 10001wall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] expand 1.5ball generate radius @rdmin number @rdmin_num box [-width*0.5] [width*0.5] ...                                                        [-height*0.5] [height*0.5]ball generate radius @rdmax number @rdmax_num box [-width*0.5] [width*0.5] ...                                                        [-height*0.5] [height*0.5]cmat default model linear method deformability emod 10e9 kratio 1.5 ball attribute density 2.7e3 damp 0.7cycle 2000 calm 50solve save sample_generate_1_1

    

最后会出现warning:

    可以看出第一个粒径会生成7004个颗粒，但是指定数目为7460，少了一点。但是当第一种颗粒占了很大区域之后，第二种颗粒目标数目为3315，却只生成了96个。如图

    针对这个问题，出现了两种解决办法——粒径膨胀法和压缩法。原理很简单，既然指定区域生成不了指定数目的颗粒，那我就减小颗粒粒径，或者增大区域。

（1） 减小颗粒粒径的方法——粒径膨胀法

    这个思路就是先将粒径减小，生成指定数目后，再将粒径增大。这里给出demo，这里定义了膨胀系数为5.0，就是先将粒径减小五倍，之后再放大。

new def par    width=1.5    height= width       rdmax=0.009    rdmin=0.006        poro=0.25     rdmax_num=math.floor((1-poro)*0.5*width*height/(math.pi*rdmax*rdmax))    rdmin_num=math.floor((1-poro)*0.5*width*height/(math.pi*rdmin*rdmin))    pengzhangxishu=5.0end@pardomain extent [-height*5] [height*5]set random 10001wall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] expand 1.5ball generate radius [rdmin/pengzhangxishu] number @rdmin_num box [-width*0.5] [width*0.5] ...                                                        [-height*0.5] [height*0.5]ball generate radius [rdmax/pengzhangxishu] number @rdmax_num box [-width*0.5] [width*0.5] ...                                                        [-height*0.5] [height*0.5]cmat default model linear method deformability emod 10e9 kratio 1.5ball attribute density 2.7e3 damp 0.7cycle 2000 calm 50cycle 4000save init_generate_pengzhangball attribute radius multiply @pengzhangxishucycle 2000 calm 50solvesave sample_generate_1_2

从console里面可以看到指定数目的颗粒生成了

模型图可以看出来颗粒非常小

之后放大粒径后的状态为：

一共10775个颗粒全部生成。

以上为粒径膨胀法的概念

（2） 增加区域面积的方法——压缩法

      压缩法的概念是，保持粒径不变，我先将区域放大，然后移动墙体，使其移动到指定的区域面积。以下为这个逻辑的demo，区域放大系数为5.0，也就是先将区域放大五倍，然后移动上部墙体，移动4.0倍的区域距离。

new def par    width=1.5    height= width       rdmax=0.009    rdmin=0.006        poro=0.25     rdmax_num=math.floor((1-poro)*0.5*width*height/(math.pi*rdmax*rdmax))    rdmin_num=math.floor((1-poro)*0.5*width*height/(math.pi*rdmin*rdmin))    pengzhangxishu=5.0    y_vel=100.0end@pardomain extent [-height*pengzhangxishu*2] [height*pengzhangxishu*2]set random 10001wall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*(pengzhangxishu-1)+height*0.5] expand 1.5ball generate radius [rdmin] number @rdmin_num box [-width*0.5] [width*0.5] ...                                            [-height*0.5] [height*(pengzhangxishu-1)+height*0.5]ball generate radius [rdmax] number @rdmax_num box [-width*0.5] [width*0.5] ...                                             [-height*0.5] [height*(pengzhangxishu-1)+height*0.5]cmat default model linear method deformability emod 10e9 kratio 1.5 ball attribute density 2.7e3 damp 0.7cycle 2000 calm 50cycle 4000save init_generate_quyuwall attribute yvel [-y_vel] range id 3solve time [height*(pengzhangxishu-1)/y_vel]wall attribute yvel 0 range id 3solvesave sample_generate_1_3

    这里先给出扩大区域内试样的状态，这里的粒径还是指定的，只是因为区域比较大，所以视觉上看颗粒显得小点。

压缩后的状态如图：

    

    可以看出来，虽然生成了指定数目的颗粒了，但是试样在上方比较密，在下方比较松，这样必然会导致试样的不均匀性，于是很多学者对其进行了拓展，基本的拓展方法为分层压缩法。

（3） 压缩法的拓展——分层压缩法

    分层压缩法的基本思路为，将压缩法分为几次（5-10次）进行压缩，这样前述的压缩法的缺点就会被忽略掉。下方给出分层压缩法的demo:

    这里用loop循环执行分层操作，每次先删除之前的顶部墙，生成新的顶部墙，然后按照压缩法进行成样，重复层数就是我们的分层压缩法。

new def par    width=1.5    height= width       rdmax=0.009    rdmin=0.006        poro=0.25         rdmax_num=math.floor((1-poro)*0.5*width*height/(math.pi*rdmax*rdmax))    rdmin_num=math.floor((1-poro)*0.5*width*height/(math.pi*rdmin*rdmin))        pengzhangxishu=5.0        y_vel=100.0        cengshu=8.0end@pardomain extent [-height*pengzhangxishu*2] [height*pengzhangxishu*2]set random 10001wall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*(pengzhangxishu-1)-height*0.5] expand 1.5cmat default model linear method deformability emod 10e9 kratio 1.5 def create_sample    loop n(1,cengshu)        dibuzuobiao=-width*0.5+height*(n-1)/cengshu        dingbuzuobiao=-width*0.5+height*(n)/cengshu                pengzhangGaoDu=dingbuzuobiao+height*(pengzhangxishu-1)/cengshu            group_string=string.build("keli_%1",n)        jieguo_string=string.build("bushu_%1",n)        command            wall delete walls range id 3            wall create id 3 vertices [-width*0.5*1.5] [pengzhangGaoDu] [width*0.5*1.5] [pengzhangGaoDu]             ball generate radius [rdmin] number [math.floor(rdmin_num/cengshu)] box [-width*0.5] [width*0.5] ...                                            [dibuzuobiao] [pengzhangGaoDu] group @group_string            ball generate radius [rdmax] number [math.floor(rdmax_num/cengshu)] box [-width*0.5] [width*0.5] ...                                            [dibuzuobiao] [pengzhangGaoDu] group @group_string            ball attribute density 2.7e3 damp 0.7            wall attribute yvel [-y_vel] range id 3            solve time [(pengzhangGaoDu-dingbuzuobiao)/y_vel]                        wall attribute yvel 0 range id 3            cycle 2000            save @jieguo_string        endcommand            endloopend@create_samplesolvesave sample_generate_1_4

这里分8次成样，下面给出执行中的试样状态：

第一次压缩：

第四次压缩：

第八次压缩：

最后状态的试样为：

这里可以看出压缩法的缺陷已经得到明显的改善。

增加压缩次数试样会更加均匀，下面为分16次压缩：

    很多学者对于压缩时候的能量传递进行研究，认为由于压缩能的影响，下部的颗粒在成样的时候应该松一点，上部的颗粒应该密一点。这时候我们可以改变孔隙率实现这个想法，具体为：如分5层，则5层的孔隙率为 1.08*poro，1.04*poro，poro，0.96*poro，0.92*poro。这个可以根据试样的情况去调。也有一些学者对孔隙率的值进行了定值研究，国内比较著名的就是分层欠压法（蒋），这部分同学们可以参考文献自己去定义，虽然我懂这个东西，但是不知道能不能讲，防止被告，还是不讲了。

二、distribute方法

（1） 基本的distribute法

    这个方法应该是generate方法后面出来的，为了克服generate方法的缺点，distribute允许颗粒有很大的重叠，这样在指定的区域内必然可以生成指定数目的颗粒。这个逻辑为：根据区域面积和孔隙率计算颗粒体积，然后根据颗粒体积计算颗粒数目，在指定区域内堆叠生成颗粒。

    为了克服大重叠引起的颗粒很大的初速度，我们需要先利用calm命令清除颗粒的动能。这里给出distribute命令的demo，这里需要计算两个粒径颗粒的孔隙率。逻辑为：Vs=(1-poro)*V*0.5，则Vsporo=(V-(1-poro)*V*0.5)/V=1-(1-poro)*0.5。

    cycle 2000 calm 50 就是为了清除颗粒的初速度，如果还是有颗粒飞出来的话，就需要使用calm 10 或者calm 5。

new def par    width=1.5    height= width       rdmax=0.009    rdmin=0.006        poro=0.25    vsporo=1-(1-poro)*0.5end@pardomain extent [-height*5] [height*5]set random 10001wall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] expand 1.5ball distribute porosity  @vsporo radius @rdmin box [-width*0.5] [width*0.5] ...                                                        [-height*0.5] [height*0.5]ball distribute porosity @vsporo radius @rdmax box [-width*0.5] [width*0.5] ...                                                        [-height*0.5] [height*0.5]cmat default model linear method deformability emod 10e9 kratio 1.5ball attribute density 2.7e3 damp 0.7cycle 2000 calm 50solve save sample_distribute_2_1

刚生成颗粒时候的状态如下，可以看出很大的重叠量

进行平衡后的试样如图：

（2） 分层成样法

    这里也是用distribute命令进行成样，只是参考分层压缩法进行操作，这样可以生成更加均匀的试样，demo如下：

new def par    width=1.5    height= width       rdmax=0.009    rdmin=0.006        poro=0.25    vsporo=1-(1-poro)*0.5    cengshu=8.0end@pardomain extent [-height*2] [height*2]set random 10001wall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] expand 1.5cmat default model linear method deformability emod 10e9 kratio 1.5def create_sample    loop n(1,cengshu)        dibuzuobiao=-width*0.5+height*(n-1)/cengshu        dingbuzuobiao=-width*0.5+height*(n)/cengshu          group_string=string.build("keli_%1",n)        jieguo_string=string.build("bushu_dist_%1",n)        command            wall delete walls range id 3            wall create id 3 vertices [-width*0.5*1.5] [dingbuzuobiao] [width*0.5*1.5] [dingbuzuobiao]            ball distribute porosity  @vsporo radius @rdmin box [-width*0.5] [width*0.5] ...                                                        [dibuzuobiao] [dingbuzuobiao] group @group_string            ball distribute porosity @vsporo radius @rdmax box [-width*0.5] [width*0.5] ...                                                        [dibuzuobiao] [dingbuzuobiao] group @group_string            ball attribute density 2.7e3 damp 0.7            cycle 2000 calm 50                        cycle 4000            save @jieguo_string        endcommand            endloopend@create_samplesolvesave sample_distribute_2_2

这里分8次distribute，每次状态为：

第一次：

第四次：

第八次：

最后状态为：

可以看出与单次distribute已经有了改善，这里分16次再算一次：

三、超级快速成样方法——Brick方法

    分层可以使得试样均匀，这个道理PFC的开发人员当然也注意到了，于是开发了Brick算法。就像其意思一样，Brick方法是先生成满足条件的一个小块，然后用这个小块拼凑起整个模型。

这个东西牛逼的是计算速度，上面的分层法计算一次基本上都需要十几分钟左右，但是brick方法只需要1秒！

Brick分为两部分：

（1）首先生成基本的砖头，注意这里的domain必须是周期边界，并且大小就是砖头的大小

（2）之后利用brick关键词进行组装

new def par    width=1.5    height= width       rdmax=0.009    rdmin=0.006        poro=0.25    vsporo=1-(1-poro)*0.5    lieshu=8.0end@pardomain extent [-width*0.5/lieshu] [width*0.5/lieshu] [-height*0.5/lieshu] [height*0.5/lieshu] condition periodicset random 10001ball distribute porosity  @vsporo radius @rdmin box [-width*0.5/lieshu] [width*0.5/lieshu] ...                                                        [-height*0.5/lieshu] [height*0.5/lieshu]ball distribute porosity @vsporo radius @rdmax box [-width*0.5/lieshu] [width*0.5/lieshu] ...                                                        [-height*0.5/lieshu] [height*0.5/lieshu]cmat default model linear method deformability emod 10e9 kratio 1.5ball attribute density 2.7e3 damp 0.7cycle 2000 calm 50solvesave brick_samplebrick make id 1ball deletedomain extent [-height*5] [height*5] condition stopwall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] expand 1.5brick assemble id 1 origin [-width*0.5] [-height*0.5] ...            size [int(lieshu)] [int(lieshu)]solvesave sample_brick

运行后的试样为：

    这里因为一个brick的数目比较少，所以效果可能并不是很理想，但是已经比直接生成效果好很多了。对于大数目的均匀试样，brick无疑是个好方法

下面给个总结吧：

总的来说，对于颗粒不是特别多的试样，分层压缩法会是一个比较好的选择；

对于颗粒很多的试样，Brick应该是第一选择