柔性簇Cluster模擬基本框架介紹

lobby

2021年10月10日 21:00

瀏覽：4897 評論：4 收藏：14

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