基于PFC與FlAC耦合的柔性三軸實(shí)驗(yàn)

lobby

2022年4月15日 14:19

瀏覽：6479 評(píng)論：14 收藏：33

前言

最近熟悉6.0最大的感受就是，之前的連續(xù)體和離散體混合模型都可以使用6.0了，這里給大家介紹一下在PFC5.0版本實(shí)現(xiàn)起來(lái)特別費(fèi)力并且效果一般的柔性三軸實(shí)驗(yàn)。

在6.0中實(shí)現(xiàn)柔性三軸方便了不是一星半點(diǎn)，而且解決了很多顆粒膜存在的問(wèn)題。有點(diǎn)實(shí)話就是：

耦合的方法使得之前的顆粒膜方法成為了笑話。

我個(gè)人用顆粒膜方法，存在的一個(gè)最大的問(wèn)題就是，顆粒膜會(huì)出現(xiàn)大的變形，因?yàn)槟さ膹椥阅Ａ勘容^小，導(dǎo)致在加載過(guò)程中顆粒膜的變形比較大。

這個(gè)對(duì)于有限元來(lái)說(shuō)問(wèn)題不大，變形大我就變大單元好了，但是離散元中的變形是通過(guò)顆粒重疊實(shí)現(xiàn)的。大變形下會(huì)使得顆粒膜顆粒之間的間距變大。這樣就會(huì)使得散體材料中的顆粒會(huì)在膜顆粒間發(fā)生逃逸。

在6.0中提供了耦合的辦法，可以使得flac中的結(jié)構(gòu)體成為pfc中的wall，這樣就可以完美解決上述問(wèn)題。

本文也參考了手冊(cè)中自帶的柔性三軸案例，在原有的單元實(shí)驗(yàn)的結(jié)構(gòu)上進(jìn)行開發(fā)。

1 成樣、預(yù)壓、圍壓

這里直接跳過(guò)。

model new def chicun_par    sample_rad=1    sample_hight=sample_rad*4        keli_rdmin=0.04    keli_rdmax=0.06end@chicun_par

model domain extent [-sample_rad*2] [sample_rad*2] [-sample_rad*2] [sample_rad*2]   [-sample_hight*1.5] [sample_hight*1.5]

[n=1.4]model random 10001wall generate cylinder base 0 0 [-sample_hight*0.5*n] axis 0 0 1 ...                    height [sample_hight*n] radius [sample_rad] resolution 0.3 cap false false

wall generate plane position 0 0 [sample_hight*0.5] dip 0 dip-direction  0wall generate plane position 0 0 [-sample_hight*0.5] dip 0 dip-direction  0

ball distribute group "shiyang" radius [keli_rdmin] [keli_rdmax] porosity 0.28 ...    range cylinder end-1 0 0 [sample_hight*0.5-keli_rdmin] ...    end-2 0 0 [-sample_hight*0.5+keli_rdmin] radius [sample_rad-keli_rdmin]

cmat default type ball-facet model linear method deform emod 100e6 kratio 1.5 property  fric 0.2cmat default type ball-ball model linear method deform emod 100e6 kratio 1.5  ball attribute density 2.7e3 damp 0.7model cycle 2000 calm 50

model solve

model save "sample"

model restore "sample"

ball property "fric" 0.5def wp_wall    wp_up=wall.find(2)    wp_down=wall.find(3)    wp_rr=wall.find(1)        loop foreach vt wall.vertexlist(wp_rr)        vert_in_ce=vt             endloopend@wp_wall

[tzz=-1e4][trr=-1e4][sevro_fac=0.5]

[do_zservo=true][do_rservo=true]



[sevro_freq=100][timestepNow=global.step-1]def servo_walls    computer_wallStress    if timestepNow<global.step then        get_gain(sevro_fac)        timestepNow+=sevro_freq    endif    if do_zservo=true then        z_vel=gz*(wszz-tzz)        wall.vel.z(wp_up)=-z_vel        wall.vel.z(wp_down)=z_vel    endif        if do_rservo=true then        r_vel_mag=(-1)*gr*(wsrr-tzz)        loop foreach vt wall.vertexlist(wp_rr)            mag=math.sqrt(wall.vertex.pos.x(vt)^2+wall.vertex.pos.y(vt)^2)            fang_normal_x=wall.vertex.pos.x(vt)/mag            fang_normal_y=wall.vertex.pos.y(vt)/mag                        r_vel=vector(fang_normal_x,fang_normal_y,0)*r_vel_mag                        wall.vertex.vel(vt)=r_vel         endloop    endif        end



def computer_chicun    x_pos=wall.vertex.pos.x(vert_in_ce)    y_pos=wall.vertex.pos.y(vert_in_ce)    wlr=math.sqrt(x_pos^2+y_pos^2)    wlz=wall.pos.z(wp_up)-wall.pos.z(wp_down)end

def computer_wallStress    computer_chicun    ding_yuanmianji=math.pi*wlr^2    wszz=(wall.force.contact.z(wp_down)-wall.force.contact.z(wp_up))*0.5/ding_yuanmianji    ce_mianji=2*math.pi*wlr*wlz    wsrr=0    loop foreach ft wall.facetlist(wp_rr)        ft_fangxiang=wall.facet.normal(ft)        loop foreach ct wall.facet.contactmap(ft)            force_in_facet=contact.force.global(ct)            wsrr+=-(math.dot(force_in_facet,ft_fangxiang))/ce_mianji               endloop    endloop    end

def get_gain(fac)    gz=0    gr=0    zonggangZ=0    zonggangR=0    loop foreach ct wall.contactmap(wp_up)        zonggangZ+=contact.prop(ct,"kn")    endloop        loop foreach ct wall.contactmap(wp_down)        zonggangZ+=contact.prop(ct,"kn")    endloop        loop foreach ct wall.contactmap(wp_rr)        zonggangR+=contact.prop(ct,"kn")    endloop        gz=fac*ding_yuanmianji/(zonggangZ*global.timestep)    gr=fac*ce_mianji/(zonggangR*global.timestep)    end





fish callback add @servo_walls -1.0

history id 1 @wszzhistory id 2 @wsrr

model cycle 1model solvemodel save "yuya"

model restore "yuya"

[tzz=-3e5][trr=-3e5]



model cycle 1 model solve

model save "weiya"

2 加柔性膜

以下代碼就可以直接取代之前柔性膜的效果，可以看出代碼量就減少了很多。

這里使用的是一個(gè)拉伸的方式建立一個(gè)圓柱體，geometry edge 通過(guò)四個(gè)1/4圓構(gòu)造一個(gè)圓，然后通過(guò)extrude 方法講圓拉伸為一個(gè)圓柱。里面的segment變量指定了節(jié)點(diǎn)的密度。

這里因?yàn)槟け容^薄，所以用的是有限元中的shell單元，shell單元有法向和切向變形模量需要指定，這里指定了切向的模量為1MPa，法向不發(fā)生變形。

之后和之前一樣，將膜分成上、中、下三部分。上、下的節(jié)點(diǎn)我們需要將速度固定住。中部的shell我們可以通過(guò)structure shell apply指定法向的應(yīng)力，這里完全替代了顆粒膜復(fù)雜繁瑣的節(jié)點(diǎn)力算法。上下的墻體我們用還是用伺服，這里沒(méi)有像傳統(tǒng)的一樣指定vmax，好像效果好很多。

這里需要注意的是，當(dāng)顆粒大小和面片大小比值超過(guò)一定量時(shí)，會(huì)使得平衡變得困難，所以也沒(méi)必要平衡到-5次方，可以根據(jù)應(yīng)力曲線平穩(wěn)后停止即可。因?yàn)榧尤嵝阅で笆綐觾?nèi)部是平衡的。

model restore "weiya"

fish callback remove @servo_walls -1.0[segments = 12][rad=wlr*1.01]

geometry edge create by-arc origin (0,0,[-wlz*0.6]) ...   start ([rad*(-1)],0,[-wlz*0.6]) end (0,[rad*(-1)],[-wlz*0.6]) ...   segments [segments]geometry edge create by-arc origin (0,0,[-wlz*0.6]) ...   start (0,[rad*(-1)],[-wlz*0.6]) end ([rad],0,[-wlz*0.6]) ...   segments [segments]geometry edge create by-arc origin (0,0,[-wlz*0.6]) ...   start ([rad],0,[-wlz*0.6]) end (0,[rad],[-wlz*0.6]) ...   segments [segments]geometry edge create by-arc origin (0,0,[-wlz*0.6]) ...   start (0,[rad],[-wlz*0.6]) end ([rad*(-1)],0,[-wlz*0.6]) ...   segments [segments];extrude the edges to make a cylindergeometry generate from-edges extrude (0,0,[wlz*1.2]) segments [segments*2]

structure shell import from-geometry 'Default' element-type dkt-cst structure shell property isotropic (1e6, 0.0) thick 0.25 density 930.0

structure node group 'top'  range position-z [wlz*0.48] [wlz*0.6]structure node group 'btm'  range position-z [-wlz*0.6] [-wlz*0.48]structure node group 'mid'  range position-z [-wlz*0.48] [wlz*0.48]structure shell group 'mid' range position-z [-wlz*0.48] [wlz*0.48]wall deletewall generate id 1 cylinder base 0 0 [wlz*0.5] axis 0 0 1 radius [wlr] height [keli_rdmax*10] one-wallwall generate id 2 cylinder base 0 0 [-wlz*0.5] axis 0 0 -1 radius [wlr] height [keli_rdmax*10] one-wall

structure damp localstructure shell apply [trr] range group 'mid'structure node fix velocity rotation range group "top"structure node fix velocity rotation range group "btm"wall servo force (0,0,[-tzz*math.pi*wlr*wlr]) activate true ...                 range id 2wall servo force (0,0,[tzz*math.pi*wlr*wlr]) activate true ...                 range id 1wall-structure create

measure create id 1 position 0 0 0  radius [wlr*0.4][mp=measure.find(1)]

def wp_wall    wp_up=wall.find(1)    wp_down=wall.find(2)end@wp_walldef jiance_measure       stressXX=measure.stress.xx(mp)    stressYY=measure.stress.yy(mp)    stressZZ=measure.stress.zz(mp)        ding_yuanmianji=math.pi*wlr^2    wszz=(wall.force.contact.z(wp_down)-wall.force.contact.z(wp_up))*0.5/ding_yuanmianjiendfish callback add @jiance_measure -1history deletehistory id 10 @stressXXhistory id 11 @stressYYhistory id 12 @stressZZhistory id 1 @wszz

model cycle 200model solve ratio-average 6e-5

model save "rouxing"

加完后的模型如圖：

基于PFC與FlAC耦合的柔性三軸實(shí)驗(yàn)的圖1

3、加載

這里和之前一樣，指定速度，需要指定給墻體，也要指定給上下的shell節(jié)點(diǎn)上。

model restore "rouxing"

wall servo activate false range id 2wall servo activate false range id 1

history delete ball attribute displacement multiply 0[strainRate=5e-2]

wall attribute velocity-z [strainRate*wlz] range id 2wall attribute velocity-z [-strainRate*wlz] range id 1



structure node initialize velocity-z [strainRate*wlz] range group "btm"structure node initialize velocity-z [-strainRate*wlz] range group "top"

[Iz0=wall.pos.z(wp_up)-wall.pos.z(wp_down)-keli_rdmax*10]def jiance       stressXX=measure.stress.xx(mp)    stressYY=measure.stress.yy(mp)    stressZZ=measure.stress.zz(mp)    wlz=wall.pos.z(wp_up)-wall.pos.z(wp_down)-keli_rdmax*10    wezz=(wlz-Iz0)/Iz0    ding_yuanmianji=math.pi*wlr^2    wszz=(wall.force.contact.z(wp_down)-wall.force.contact.z(wp_up))*0.5/ding_yuanmianjiendfish callback add @jiance -1

history id 10 @stressXXhistory id 11 @stressYYhistory id 12 @stressZZhistory id 1 @wszzhistory id 2 @wezz

[stop_me=0]def stop_me    if wezz<-20e-2 then        stop_me=1    endifend

[baocunpinlv=1e-2][time_record=wezz+1][count=0]def savefile        if time_record-wezz >= baocunpinlv then        filename=string.build("jieguo_%1",count)        command           model save @filename        endcommand        time_record=wezz        count +=1    endif    endfish callback add @savefile -1.0

model solve fishhalt @stop_memodel save "result"