指定级配的分层压缩法成样

2年前浏览2936

前述文章讲述了基本的成样方法（简述PFC成样方法），其中通过对比也能感受到分层压缩法在均匀试样形成中的优越性。

前面介绍的成样方法中，试样只有两种粒径，其实已经讲述了如何根据级配计算颗粒数的方法。这里为了方便使用，给出指定级配的试样的分层压缩方法。

一、二维矩形试样

下面为参数部分：

new def par    width=0.3    height= width*2         y_vel=10.0        cengshu=8.0    pengzhangxishu=5.0    tbjipei=table.create("jipei")    table(tbjipei,0.004)=0.1    table(tbjipei,0.005)=0.12    table(tbjipei,0.0055)=0.14    table(tbjipei,0.006)=0.3    table(tbjipei,0.0065)=0.6    table(tbjipei,0.007)=0.8    table(tbjipei,0.0075)=0.95    table(tbjipei,0.008)=1     poro=0.12    ;dengchaxian=0.02end@par

这里使用table记录级配，x和y分别是粒径和对应粒径的质量百分比。注意这里的粒径没有连续的变化，都是突变的，如果粒径分的足够细的话也就是曲线了。基本的参数还有尺寸以及层数，还有每次压缩时候区域的膨胀系数。

table可以显示出级配曲线：

之后进行分层压缩，这里预留了每层孔隙率变化的代码部分，读者可以自己体会一下注释部分。



new def par    width=0.3    height= width*2         y_vel=10.0        cengshu=8.0    pengzhangxishu=5.0    tbjipei=table.create("jipei")    table(tbjipei,0.004)=0.1    table(tbjipei,0.005)=0.12    table(tbjipei,0.0055)=0.14    table(tbjipei,0.006)=0.3    table(tbjipei,0.0065)=0.6    table(tbjipei,0.007)=0.8    table(tbjipei,0.0075)=0.95    table(tbjipei,0.008)=1     poro=0.12    ;dengchaxian=0.02end@pardomain extent [-height*10] [height*10]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 100e6 kratio 1.5 

def create_sample    tb_ceng=table.create("ceng_poro")    loop n(1,cengshu)        dibuzuobiao=-height*0.5 height*(n-1)/cengshu        dingbuzuobiao=-height*0.5 height*(n)/cengshu                pengzhangGaoDu=dingbuzuobiao height*(pengzhangxishu-1)/cengshu           command            wall delete walls range id 3            wall create id 3 vertices [-width*0.5*1.5] [pengzhangGaoDu] [width*0.5*1.5] [pengzhangGaoDu]         endcommand        poroceng=poro;-(n-(cengshu 1)*0.5)*dengchaxian/(cengshu-1)        table(tb_ceng,n)=poroceng          keligeoup=table.size(tbjipei)        Vs=height*width*(1-poroceng)/cengshu         cenggName=string.build("ceng%1",n)        loop m(1,keligeoup)            keligid=keligeoup 1-m            vslocal=0            if keligid=1 then                vslocal=Vs*table.y(tbjipei,keligid)            else                            vslocal=Vs*(table.y(tbjipei,keligid)-table.y(tbjipei,keligid-1))            endif            kelinum=vslocal/(math.pi*(table.x(tbjipei,keligid)*0.5)^2)            command                 ball generate radius [table.x(tbjipei,keligid)*0.5] number [math.ceiling(kelinum)] box [-width*0.5] [width*0.5] ...                    [dibuzuobiao] [pengzhangGaoDu] group @cenggName tries 20000                ball attribute density 2.7e3 damp 0.7            endcommand        endloop                 command            wall attribute yvel [-y_vel] range id 3            solve time [(pengzhangGaoDu-dingbuzuobiao)/y_vel]                        wall attribute yvel 0 range id 3            cycle 2000        endcommand            endloopend@create_samplesolvesave sample

成好的试样为：

这里分了八层成样，对试样采用测量圆检测其孔隙率，方法在前述的成样方法中已讲述（简述PFC成样方法）

测量圆布置为：

640 (1).png

孔隙率分布为：

可以看出孔隙率与预计的要差一点，这是因为计算粒径的时候没有考虑到重叠，所以测量圆测的孔隙率要小一点。这个在成密样的时候会出现这个问题，所以可能之后的孔隙率要预设大一点。

二、三维方形试样

这里道理是差不多的，只是多了一个维度而已：



new def par    width=0.1    height=width*2.0         y_vel=10.0        cengshu=8.0    pengzhangxishu=5.0    tbjipei=table.create("jipei")    table(tbjipei,0.004)=0.1    table(tbjipei,0.005)=0.12    table(tbjipei,0.0055)=0.14    table(tbjipei,0.006)=0.3    table(tbjipei,0.0065)=0.6    table(tbjipei,0.007)=0.8    table(tbjipei,0.0075)=0.95    table(tbjipei,0.008)=1         poro=0.25    ;dengchaxian=0.02end@pardomain extent [-height*pengzhangxishu*2] [height*pengzhangxishu*2]set random 10001wall generate box [-width*0.5] [width*0.5] [-width*0.5] [width*0.5] [-height*0.5] [height*(pengzhangxishu-1)-height*0.5] expand 1.5cmat default model linear method deformability emod 100e6 kratio 1.5 

def create_sample    tb_ceng=table.create("ceng_poro")    loop n(1,cengshu)        dibuzuobiao=-height*0.5 height*(n-1)/cengshu        dingbuzuobiao=-height*0.5 height*(n)/cengshu                pengzhangGaoDu=dingbuzuobiao height*(pengzhangxishu-1)/cengshu           command            wall delete walls range id 2            wall create id 2 vertices [-width*0.5*1.5] [-width*0.5*1.5] [pengzhangGaoDu] ...                                       [width*0.5*1.5] [width*0.5*1.5] [pengzhangGaoDu] ...                                      [width*0.5*1.5] [-width*0.5*1.5] [pengzhangGaoDu] ...                                       [-width*0.5*1.5] [-width*0.5*1.5] [pengzhangGaoDu] ...                                       [width*0.5*1.5] [width*0.5*1.5] [pengzhangGaoDu] ...                                      [-width*0.5*1.5] [width*0.5*1.5] [pengzhangGaoDu]        endcommand        poroceng=poro;-(n-(cengshu 1)*0.5)*dengchaxian/(cengshu-1)        table(tb_ceng,n)=poroceng          keligeoup=table.size(tbjipei)        Vs=height*width*width*(1-poroceng)/cengshu         cenggName=string.build("ceng%1",n)        loop m(1,keligeoup)            keligid=keligeoup 1-m            vslocal=0            if keligid=1 then                vslocal=Vs*table.y(tbjipei,keligid)            else                            vslocal=Vs*(table.y(tbjipei,keligid)-table.y(tbjipei,keligid-1))            endif            kelinum=vslocal/((4.0/3.0)*math.pi*(table.x(tbjipei,keligid)*0.5)^3)            command                 ball generate radius [table.x(tbjipei,keligid)*0.5] number [math.ceiling(kelinum)] box [-width*0.5] [width*0.5] [-width*0.5] [width*0.5]  ...                    [dibuzuobiao] [pengzhangGaoDu] group @cenggName tries 2000000                 ball attribute density 2.7e3 damp 0.7            endcommand        endloop                 command            wall attribute zvel [-y_vel] range id 2            solve time [(pengzhangGaoDu-dingbuzuobiao)/y_vel]                        wall attribute zvel 0 range id 2            cycle 2000        endcommand            endloopend@create_samplesolveball delete range x [-width*0.5] [width*0.5] not y [-width*0.5] [width*0.5] not z [-height*0.5] [height*0.5] notsave sample_lifang

这里为了方便计算，试样级配稍微紧一点：

最后结果为：

测量圆布置为：

检测结果为：

和二维一样，孔隙率会略低。

三、三维圆柱试样

这里算的时候有两个问题：

1）圆柱形内使用generate时，数目一直不够，于是我先方形generate，然后删除圆柱其余部分。不知道是不是版本问题，直接圆柱形generate生成的代码也保留在煮熟中，高版本的读者可以尝试一下。

2）圆柱体无摩擦平衡一直平衡不下来，经过检查，可能是可以facet很多有关系，将圆柱墙的facet减少，会容易平衡一点，但是圆柱形状也出不来，现在我也没想好比较好的办法，只能成样的时候加了fric加快平衡。



new def par    radshiyang=0.04    height= radshiyang*4        y_vel=100.0        cengshu=8.0    pengzhangxishu=5.0    tbjipei=table.create("jipei")    table(tbjipei,0.002)=0.1    table(tbjipei,0.003)=0.14    table(tbjipei,0.004)=0.3    table(tbjipei,0.005)=0.6    table(tbjipei,0.006)=0.8    table(tbjipei,0.009)=1     poro=0.25    dengchaxian=0.02end@pardomain extent [-height*pengzhangxishu*2] [height*pengzhangxishu*2]set random 10001wall generate id 3 cylinder base 0 0 [-height*0.5*pengzhangxishu] axis 0 0  1 ...        height [height*pengzhangxishu] radius [radshiyang]  cap false false wall create id 1 vertices [-radshiyang*1.5] [-radshiyang*1.5] [-height*0.5] ...                                       [radshiyang*1.5] [radshiyang*1.5] [-height*0.5] ...                                      [radshiyang*1.5] [-radshiyang*1.5] [-height*0.5] ...                                       [-radshiyang*1.5] [-radshiyang*1.5] [-height*0.5] ...                                       [radshiyang*1.5] [radshiyang*1.5] [-height*0.5] ...                                      [-radshiyang*1.5] [radshiyang*1.5] [-height*0.5] cmat default model linear method deformability emod 10e9 kratio 1.5 

def create_sample    tb_ceng=table.create("ceng_poro")    loop n(1,cengshu)        dibuzuobiao=-height*0.5 height*(n-1)/cengshu        dingbuzuobiao=-height*0.5 height*(n)/cengshu                pengzhangGaoDu=dingbuzuobiao height*(pengzhangxishu-1)/cengshu           command            wall delete walls range id 2            wall create id 2 vertices [-radshiyang*1.5] [-radshiyang*1.5] [pengzhangGaoDu] ...                                       [radshiyang*1.5] [radshiyang*1.5] [pengzhangGaoDu] ...                                      [radshiyang*1.5] [-radshiyang*1.5] [pengzhangGaoDu] ...                                       [-radshiyang*1.5] [-radshiyang*1.5] [pengzhangGaoDu] ...                                       [radshiyang*1.5] [radshiyang*1.5] [pengzhangGaoDu] ...                                      [-radshiyang*1.5] [radshiyang*1.5] [pengzhangGaoDu]        endcommand        poroceng=poro-(n-(cengshu 1)*0.5)*dengchaxian/(cengshu-1)        table(tb_ceng,n)=poroceng          keligeoup=table.size(tbjipei)        Vs=(math.pi*radshiyang^2)*height*(1-poroceng)/cengshu         cenggName=string.build("ceng%1",n)        loop m(1,keligeoup)            keligid=keligeoup 1-m            vslocal=0            if keligid=1 then                vslocal=Vs*table.y(tbjipei,keligid)            else                            vslocal=Vs*(table.y(tbjipei,keligid)-table.y(tbjipei,keligid-1))            endif            ;kelinum=vslocal/((4.0/3.0)*math.pi*(table.x(tbjipei,keligid)*0.5)^3)            kelinum=vslocal*4/math.pi/((4.0/3.0)*math.pi*(table.x(tbjipei,keligid)*0.5)^3)            command                ball generate radius [table.x(tbjipei,keligid)*0.5] number [math.ceiling(kelinum)] box [-radshiyang] [radshiyang] [-radshiyang] [radshiyang]  ...                    [dibuzuobiao] [pengzhangGaoDu] group @cenggName tries 2000000                ball delete range cylinder ...                 end1 0 0 [-height*0.5]  end2 0 0 [pengzhangGaoDu] radius [radshiyang*0.96] not                ; ball generate radius [table.x(tbjipei,keligid)*0.5] number [math.ceiling(kelinum)] ...;                 group @cenggName tries 2000000 range cylinder ...;                 end1 0 0 [dibuzuobiao]  end2 0 0 [pengzhangGaoDu] radius [radshiyang]                  ball attribute density 2.7e3 damp 0.7            endcommand        endloop                 command            wall attribute zvel [-y_vel] range id 2            solve time [(pengzhangGaoDu-dingbuzuobiao)/y_vel]                        wall attribute zvel 0 range id 2            cycle 2000        endcommand            endloopend@create_sampleball delete range cylinder ...                 end1 0 0 [-height*0.5]  end2 0 0 [height*0.5] radius [radshiyang] notcmat default type ball-ball model linear method deformability emod 10e9 kratio 1.5 property fric 0.2cmat default type ball-facet model linear method deformability emod 10e9 kratio 1.5 property fric 0.2cmat applysolve



save sample_yuanxing