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

2年前浏览6549

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

第二步是进行预压：

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的时候加。这一步主要是对弹性参数进行调整。

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

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

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

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

二、进行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

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

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

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