柔性簇Cluster模拟基本框架介绍

2021年10月10日 浏览：2852 评论：3 收藏：8

    

    PFC中基本单元有ball和clump，当然最新的有block。ball模拟圆形刚性颗粒，clump模拟不规则形状刚性颗粒。对于非刚性颗粒来说，柔性簇Cluster方法便可以比较好的去模拟了。但是目前来说cluster方法介绍的并不是很多，这节课主要针对于cluster的模拟框架进行介绍。

    本文介绍的cluster框架有两个主要步骤：1、导出cluster中ball的位置信息；2、在双轴或者别的实验的框架中将ball颗粒换为cluster。

    本文只介绍圆形柔性簇的方法，非圆形道理是一样的。

一、标定cluster单元的微观参数

    首先在按形状生成墙，并且在形状内生成颗粒。这里注意形状中心放在原点，这样可以在预压的时候进行径向伺服。

newset random 10001def par    lijing=0.1        rdmax=0.009    rdmin=0.006        poro=0.08end@pardomain extent [-lijing*2] [lijing*2]wall generate circle position 0 0  radius [lijing*0.5] ball distribute radius @rdmin @rdmax porosity @poro ...             range annulus center 0 0 radius 0 [lijing*0.5] ball attribute density 2.7e3 damp 0.7
cmat default model linear method deform emod 1e8 kratio 1.5
cycle 1000 calm 2cmat default model linear method deform emod 1e8 kratio 1.5 property fric 0.2cmat applysolvesave init_model

柔性簇Cluster模拟基本框架介绍的图1

   

 第二步是进行预压：

restore init_model[trr=-1e6][sevro_factor=0.2]
[timestep_now=global.step-1][gr_freq=100]def sevro_walls    computer_stress    if global.step>timestep_now then        get_gr(sevro_factor)        timestep_now+=gr_freq    endif    rVel=gr*(wrss-trr) ;wrss小于trr，rvel是正值，墙收缩，速度和方向矢量反向    loop foreach vt wall.vertex.list         vt_pos=wall.vertex.pos(vt)         fangxiang_mag=math.sqrt(comp.x(vt_pos)^2+Comp.y(vt_pos)^2)                fang_normal=vector(comp.x(vt_pos),comp.y(vt_pos))/fangxiang_mag        vel=rVel*fang_normal*(-1)        wall.vertex.vel(vt) = vel     endloop    end[wp=wall.find(1)][vt1=wall.vertex.find(1) ]def get_chicun    pos=wall.vertex.pos(vt1)     wlr=math.sqrt(comp.x(pos)^2+Comp.y(pos)^2)end
def computer_stress    get_chicun    wrss=0    loop foreach ct wall.contactmap(wp)        ct_force=contact.force.global(ct)        force_mag=math.sqrt(comp.x(ct_force)^2+Comp.y(ct_force)^2)        wrss+=force_mag    endloop        wrss=-wrss/(2*math.pi*wlr)    end
def get_gr(fac)    zonggangR=0    loop foreach ct wall.contactmap(wp)        zonggangR+=contact.prop(ct,"kn")    endloop    gr=fac*2*math.pi*wlr/(zonggangR*global.timestep)end@get_gr(@sevro_factor)
set fish callback -1.0 @sevro_walls
history id 1 @wrsshistory id 2 @grcycle 1solve aratio 1e-5
save yuya

这个概念和岩石差不多，调整到合适的赋存应力。

柔性簇Cluster模拟基本框架介绍的图2

第三步是给参数了，这里的参数用的常用的煤岩的参数。这里注意不需要给胶结，我们需要的是加胶结前颗粒的位置，胶结可以在模拟cluster的时候加。这一步主要是对弹性参数进行调整。

restore yuya[pb_coh=10e6][ten_coh=2.7]
cmat default type ball-facet model linear method deformability emod 10e8 kratio 1.5 property fric 0.2

cmat default type ball-ball model linearpbond method deformability ...            emod 12e8 kratio 1.5 pb_deformability emod 54e8 kratio 1.5 ...    property pb_coh [pb_coh] pb_ten [pb_coh*ten_coh] pb_fa 50 fric 0cmat applyclean
cycle 1solvesave fucan

柔性簇Cluster模拟基本框架介绍的图3

    最后一步就是导出颗粒的位置了，这里使用table实现数据的传输，将颗粒的位置和半径存入table中。

restore fucandef create_table    tb_pos=table.create("cluster_template_pos_big")    loop foreach bp ball.list         table(tb_pos,ball.pos.x(bp))=ball.pos.y(bp)    endloop        tb_rad=table.create("cluster_template_rad_big")    loop foreach bp ball.list         table(tb_rad,ball.pos.x(bp))=ball.radius(bp)    endloopend@create_tabletable cluster_template_pos_big write cluster_template_pos_big truncate table cluster_template_rad_big write cluster_template_rad_big truncate def get_lijing    cluster_lijing=0    loop foreach bp ball.list        lijing_temp=math.mag(ball.pos(bp))+ball.radius(bp)        if lijing_temp*2.0> cluster_lijing then            cluster_lijing=lijing_temp*2.0        endif    endloopend@get_lijingsave ball_pos

这里注意一下颗粒的粒径，由于预压的时候圆形半径会发生改变，一开始是0.1m，后面的0.123才是cluster真实的粒径。这里需要注意一下。

同样的道理可以再生成一个粒径的cluster模板，目录下生成了四个table，这个是我们后面需要用到的。

柔性簇Cluster模拟基本框架介绍的图4

二、进行cluster双轴

    框架还是用我常用的双轴框架。这里注意粒径的选取取决于上一步做的cluster模板。我们这里只有两个粒径。

1、成样：

new

set random 10001

def par    width=1.5    lijing1=0.123    lijing2=0.09443    height=width*2.0        poro=0.12      end@pardomain extent [-width*5] [width*5]wall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] expand 1.5ball distribute porosity @poro numbins 2 ...    bin 1 radius [lijing1*0.5] v 0.5 group cluster_template_big ...    bin 2 radius [lijing2*0.5] v 0.5 group cluster_template_small ...    box [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] ball attribute density 2.7e3 damp 0.7cmat default model linear method deformability emod 10e8 kratio 1.5 cycle 2000 calm 5solvesave sample

柔性簇Cluster模拟基本框架介绍的图5



2、预压：

restore sample

[txx=-1e4][tyy=-1e4][sevro_factor=0.5][do_xSevro=true][do_ySevro=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)        wall.vel.x(wpRight)=-Xvel        wall.vel.x(wpLeft)=Xvel    endif        if do_ySevro=true then        Yvel=gy*(wyss-tyy)        wall.vel.y(wpUp)=-Yvel        wall.vel.y(wpDown)=Yvel    endifend

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

def computer_chiCun    wlx=wall.pos.x(wpRight)-wall.pos.x(wpLeft)    wly=wall.pos.y(wpUp)-wall.pos.y(wpDown)end

def compute_stress    computer_chiCun    wxss=-(wall.force.contact.x(wpRight)-wall.force.contact.x(wpLeft))*0.5/wly    wyss=-(wall.force.contact.y(wpUp)-wall.force.contact.y(wpDown))*0.5/wlxend

def get_g(fac)    gx=0    gy=0    zongKNX=100e6*2*10    zongKNY=100e6*2*10    loop foreach ct wall.contactmap(wpLeft)        zongKNX+=contact.prop(ct,"kn")    endloop    loop foreach ct wall.contactmap(wpRight)        zongKNX+=contact.prop(ct,"kn")    endloop    loop foreach ct wall.contactmap(wpUp)        zongKNY+=contact.prop(ct,"kn")    endloop    loop foreach ct wall.contactmap(wpDown)        zongKNY+=contact.prop(ct,"kn")    endloop    gx=fac*wly/(zongKNX*global.timestep)    gy=fac*wlx/(zongKNY*global.timestep)    end

set fish callback -1.0 @sevro_walls
history id 1 @wxsshistory id 2 @wyss
history id 3 @gxhistory id 4 @gy
cycle 1

solve save yuya

3、将ball换为cluster

    这里注意将第一步生成的table复制粘贴到这个模拟的目录下，导入后进行cluster的生成。这里面细节比较多，就不一一讲解了，读者可以自己理解一下。

restore yuya

cmat default type ball-facet model linear method deformability emod 10e8 kratio 1.5 

[pb_coh=10e6][ten_coh=2.7]contact groupbehavior andcmat default type ball-ball model linearpbond method deformability ...            emod 12e8 kratio 1.5 pb_deformability emod 54e8 kratio 1.5 ...    property pb_coh [pb_coh] pb_ten [pb_coh*ten_coh] pb_fa 50 fric 0.5 table cluster_template_pos_big read cluster_template_pos_bigtable cluster_template_rad_big read cluster_template_rad_big

table cluster_template_pos_small read cluster_template_pos_smalltable cluster_template_rad_small read cluster_template_rad_small[tb_pos_big=table.find("cluster_template_pos_big")][tb_rad_big=table.find("cluster_template_rad_big")][tb_pos_small=table.find("cluster_template_pos_small")][tb_rad_small=table.find("cluster_template_rad_small")]

[tb_size_big=table.size(tb_pos_big)][tb_size_small=table.size(tb_pos_small)]def genghuan_cluster    count_cluster=0    loop foreach bp ball.groupmap("cluster_template_big")        x_pos=ball.pos.x(bp)        y_pos=ball.pos.y(bp)        ball.delete(bp)        jiaodu=math.random.uniform*math.pi*2        sin_jiaodu=math.sin(jiaodu)        cos_jiaodu=math.cos(jiaodu)        loop n(1,tb_size_big)            ct_posX=table.x(tb_pos_big,n)            ct_posY=table.y(tb_pos_big,n)            ct_rad=table.y(tb_rad_big,n)            cluster_pos_x=x_pos+(-1)*sin_jiaodu*ct_posX+cos_jiaodu*ct_posY            cluster_pos_y=y_pos+cos_jiaodu*ct_posX+sin_jiaodu*ct_posY            pos_vec=vector(cluster_pos_x,cluster_pos_y)            bp_temp=ball.create(ct_rad,pos_vec)            ball.group(bp_temp)="big_cluster"                       geoup_2_name=string.build("cluter_num_%1",count_cluster)            ball.group(bp_temp,2)=geoup_2_name        endloop        command            clean            contact method bond gap [lijing1*0.08*0.2] range group @geoup_2_name slot 2                     endcommand        count_cluster+=1    endloop        loop foreach bp ball.groupmap("cluster_template_small")        x_pos=ball.pos.x(bp)        y_pos=ball.pos.y(bp)        ball.delete(bp)        jiaodu=math.random.uniform*math.pi*2        sin_jiaodu=math.sin(jiaodu)        cos_jiaodu=math.cos(jiaodu)        loop n(1,tb_size_small)            ct_posX=table.x(tb_pos_small,n)            ct_posY=table.y(tb_pos_small,n)            ct_rad=table.y(tb_rad_small,n)            cluster_pos_x=x_pos+(-1)*sin_jiaodu*ct_posX+cos_jiaodu*ct_posY            cluster_pos_y=y_pos+cos_jiaodu*ct_posX+sin_jiaodu*ct_posY            pos_vec=vector(cluster_pos_x,cluster_pos_y)            bp_temp=ball.create(ct_rad,pos_vec)            ball.group(bp_temp)="small_cluster"                    geoup_2_name=string.build("cluter_num_%1",count_cluster)            ball.group(bp_temp,2)=geoup_2_name        endloop        command            clean            contact method bond gap [lijing1*0.08*0.2] range group @geoup_2_name slot 2                     endcommand        count_cluster+=1    endloopend@genghuan_clusterball attribute density 2.7e3 damp 0.7

cycle 1 solvesave cluster_sample

柔性簇Cluster模拟基本框架介绍的图6



    可以看一下pb_state观察一下这其中的效果。可以看到不同cluster之间是没有胶结的。


柔性簇Cluster模拟基本框架介绍的图7



4、加围压


    后面就跟双轴一样了，只是第三步不一样。


restore cluster_sample[txx=-300e3][tyy=-300e3]

cycle 1solvesave weiya

5、加载：

restore weiyaset mech age 0ball attribute displacement multiply 0[streainRatio=1e-2]

[do_xSevro=true][do_ySevro=false]wall attribute yvelocity [wly*streainRatio*0.5] range id 1wall attribute yvelocity [-wly*streainRatio*0.5] range id 3

[Ix0=wlx][Iy0=wly]def computer_strain    wexx=(wlx-Ix0)/Ix0    weyy=(wly-Iy0)/Iy0    weVol=wexx+weyy       pianyingli=(wyss-wxss)end

set fish callback -1.1 @computer_strain

history deletehistory id 1 @wysshistory id 2 @weyyhistory id 3 @wxsshistory id 4 @wexxhistory id 5 @weVolhistory id 6 @pianyingli
[stop_me=0]def stop_me    if weyy<-20e-2 then        stop_me=1    endifend
[baocunpinlv=-1e-2][time_record=weyy+1][count=0]def savefile        if weyy-time_record <= baocunpinlv then        filename=string.build("jieguo_%1",count)        command            save @filename        endcommand        time_record=weyy        count +=1    endif    endset fish callback -1.0 @savefile
solve fishhalt @stop_mesave result

看一下应力应变，由于颗粒不多，波动还是比较大的。

柔性簇Cluster模拟基本框架介绍的图8

位移图和刚性颗粒差不多。

柔性簇Cluster模拟基本框架介绍的图9

我们看一下细观：

柔性簇Cluster模拟基本框架介绍的图10


对于低围压来说，基本上都不会坏，这里的特性可能和刚性差不多。


这里尝试一下1mpa的围压：


应力应变：


柔性簇Cluster模拟基本框架介绍的图11


细观：

可以看到很多颗粒被压坏了。

柔性簇Cluster模拟基本框架介绍的图12


放大点看看：


柔性簇Cluster模拟基本框架介绍的图13


    本文只给出模拟框架，分析方面不多，读者可以根据自己需求加入一些后处理。



(3条)
默认 最新
赞👍🏻
评论 点赞
赞
评论 点赞

查看更多评论 >

点赞 4 评论 3 收藏 8
关注