首页/文章/ 详情

基于PFC与FlAC耦合的柔性三轴实验

lobby

2年前浏览6875

前言

最近熟悉6.0最大的感受就是，之前的连续体和离散体混合模型都可以使用6.0了，这里给大家介绍一下在PFC5.0版本实现起来特别费力并且效果一般的柔性三轴实验。

在6.0中实现柔性三轴方便了不是一星半点，而且解决了很多颗粒膜存在的问题。有点实话就是：

耦合的方法使得之前的颗粒膜方法成为了笑话。

我个人用颗粒膜方法，存在的一个最大的问题就是，颗粒膜会出现大的变形，因为膜的弹性模量比较小，导致在加载过程中颗粒膜的变形比较大。

这个对于有限元来说问题不大，变形大我就变大单元好了，但是离散元中的变形是通过颗粒重叠实现的。大变形下会使得颗粒膜颗粒之间的间距变大。这样就会使得散体材料中的颗粒会在膜颗粒间发生逃逸。

在6.0中提供了耦合的办法，可以使得flac中的结构体成为pfc中的wall，这样就可以完美解决上述问题。

本文也参考了手册中自带的柔性三轴案例，在原有的单元实验的结构上进行开发。

1 成样、预压、围压

这里直接跳过。

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 加柔性膜

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

这里使用的是一个拉伸的方式建立一个圆柱体，geometry edge 通过四个1/4圆构造一个圆，然后通过extrude 方法讲圆拉伸为一个圆柱。里面的segment变量指定了节点的密度。

这里因为膜比较薄，所以用的是有限元中的shell单元，shell单元有法向和切向变形模量需要指定，这里指定了切向的模量为1MPa，法向不发生变形。

之后和之前一样，将膜分成上、中、下三部分。上、下的节点我们需要将速度固定住。中部的shell我们可以通过structure shell apply指定法向的应力，这里完全替代了颗粒膜复杂繁琐的节点力算法。上下的墙体我们用还是用伺服，这里没有像传统的一样指定vmax，好像效果好很多。

这里需要注意的是，当颗粒大小和面片大小比值超过一定量时，会使得平衡变得困难，所以也没必要平衡到-5次方，可以根据应力曲线平稳后停止即可。因为加柔性膜前式样内部是平衡的。

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"

加完后的模型如图：

3、加载

这里和之前一样，指定速度，需要指定给墙体，也要指定给上下的shell节点上。

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"