PFC砂土柔性双轴

2年前浏览4595

前言：

柔性双轴和三轴的代码一年半以前就已经写好框架了，但有一些Bug，一直没有时间去整理。最近计算力有点空闲，刚好今天有点酒意，将柔性双轴的代码先整理了一下，基本上土的力学特性是可以反应了，剩下的就靠各位去继续完善了。

首先我们得搞清楚为什么要做柔性双轴或者三轴。真三轴其实是更加符合土单元概念的力学试验，但是在现实中真三轴的难度可以说是假三轴的百倍以上。这就导致了目前很多土力学实验都是假三轴。而我们做参数标定，是需要做和现实一致的单元实验，调整微观参数使其宏观特性一致的，所以数值模拟中做假三轴比真三轴更好。

数值模拟是在还原现实的基础上去研究更多力学特性，所以第一步还原现实我们需要首先做到位。

一、成样和预压

这里还是和之前的一样，成样代码如下：

newdomain extent -5 5set random 10001def par    width=1.5    height=width*2    poro=0.18end@parwall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] expand 1.5cmat default model linear method deformability emod 100e6 kratio 1.5 ball distribute  radius 0.006 0.009 porosity @poro ...         range x [-width*0.5 0.001]  [width*0.5-0.001] ...               y [-height*0.5 0.001] [height*0.5-0.001]ball attribute density 7800 damp 0.7cycle 1000 calm 10solve ball delete range x [-width] [-width*0.5]  ball delete range x [width*0.5] [width]ball delete range y [-height] [-height*0.5] ball delete range y [height*0.5] [height]

save ball_sample

预压代码如下，注意在压之前给参数，这是我认为比较好的砂土模拟步骤。

restore ball_sample[fric_shiyang=0.5][rrfric_shiyang=0.2][emod_shiyang=100e6]ball group shiyang cmat add 1 model rrlinear method deformability emod @emod_shiyang kratio 1.5 property fric @fric_shiyang rr_fric @rrfric_shiyang  ...                                        range group  shiyangcmat apply   [txx=-2e4][tyy=-2e4]

[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;   sudu        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@compute_stress

def get_g(fac)    computer_chiCun    gx=0    gy=0    zongKNX=100e6*10    zongKNY=100e6*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

@compute_stress

set fish callback -1.0 @sevro_walls

history id 1 @wxsshistory id 2 @wysscycle 1solve aratio 1e-5save yuya

二、柔性膜

这就是本文的重点了，下图为柔性膜颗粒计算的基本逻辑，静水压力为stressW，这里将颗粒间看作是连续的膜，膜的长度就是接触的长度L，于是静水压力作用在整个膜上的力为fw。我们再将这个力分解到膜的两端，也就是接触两端的颗粒，分别为0.5*fw。

我们首先生成上下的加载板，这里用vertice生成异形墙体就行，这里采用指定坐标的方法。

之后生成膜颗粒，这里循环生成即可，上下加载板中都会有膜颗粒，以模拟膜绑在加载板上的样子。

这两步走完后是这个样子：

之后给膜颗粒加接触，这里指定模量为7e6，也是一般橡胶的模量，强度设为1e300，这样膜是不会坏的。

下面就是最重要的施加膜颗粒上的作用力了。

disp就是接触的长度L，fmag就是0.5*fw。

dirt为两个膜颗粒的方向，注意这里我们需要旋转90度，逆时针还是顺时针取决于是左边的膜还是右边的膜，以及计算时候bp2在bp的上方还是下方。

左边和右边我们在生成的时候已经分好组了。上方下方我们可以根据颗粒的id来确定。

举个例子：如果是左边膜，然后bp2在bp的上方，根据bp和bp2计算得到的力方向应当是dirt顺时针旋转90度。

至于旋转矩阵就是高中知识了，(x,y)逆时针旋转90度为(-y,x)。

[weiya=math.abs(txx)]def upadat_bianjie_force    loop foreach bp ball.groupmap("left_bianjie")        if ball.contactnum(bp)>=2 then            f=vector(0,0)            loop foreach ct ball.contactmap(bp)                if contact.model(ct)="linearcbond" then                    bp2=contact.end1(ct)                    if bp2=bp then                        bp2=contact.end2(ct)                    endif                    dirt=ball.pos(bp2)-ball.pos(bp)                    disp=math.sqrt(comp.x(dirt)^2 comp.y(dirt)^2)                    fmag=weiya*disp/2.0                     factor=1                    if ball.id(bp2)>ball.id(bp) then                        factor=-1                    endif                     norm=factor*vector(comp.y(dirt)*(-1),comp.x(dirt))/disp                     f= f norm*fmag                endif             endloop            ball.force.app(bp)=f        endif

    endloop     loop foreach bp ball.groupmap("right_bianjie")        if ball.contactnum(bp)>=2 then            f=vector(0,0)            loop foreach ct ball.contactmap(bp)                if contact.model(ct)="linearcbond" then                    bp2=contact.end1(ct)                    if bp2=bp then                        bp2=contact.end2(ct)                    endif                    dirt=ball.pos(bp2)-ball.pos(bp)                    disp=math.sqrt(comp.x(dirt)^2 comp.y(dirt)^2)                    fmag=weiya*disp/2.0                    factor=-1                    if ball.id(bp2)>ball.id(bp) then                        factor=1                    endif                    norm=factor*vector(comp.y(dirt)*(-1),comp.x(dirt))/disp                      f= f norm*fmag                endif            endloop            ball.force.app(bp)=f        endif

    endloopend

另外上下的加载板也需要伺服，我们还需要将加载板上的颗粒固定住。



def initUpDownMo     loop foreach bp ball.groupmap("left_bianjie")        if ball.pos.y(bp)>wly*0.5-bianjie_rad then            ball.group(bp,10)="upMo"         endif        if ball.pos.y(bp)<-wly*0.5 bianjie_rad then            ball.group(bp,10)="downMo"         endif     endloop     loop foreach bp ball.groupmap("right_bianjie")        if ball.pos.y(bp)>wly*0.5-bianjie_rad then            ball.group(bp,10)="upMo"        endif        if ball.pos.y(bp)<-wly*0.5 bianjie_rad then            ball.group(bp,10)="downMo"                    endif     endloopend@initUpDownMo

ball fix vel spin range group upMo slot 10ball fix vel spin range group downMo slot 10

这里的upMO和downMo就是加载板上的颗粒。

之后就是伺服的过程了，不去多讲。

完整柔性膜部分的代码如下：

restore yuya[tyy=-10e4][txx=-10e4] set fish callback -1.0 remove @sevro_walls[bianjie_rad=0.002][freq=200]wall deletewall create vertices [-wlx*0.5] [wly*0.5*1.5] [-wlx*0.5] [wly*0.5] ...                        [-wlx*0.5] [wly*0.5] [wlx*0.5] [wly*0.5] ...                        [wlx*0.5] [wly*0.5] [wlx*0.5] [wly*0.5*1.5] id 1wall create vertices [-wlx*0.5] [-wly*0.5*1.5] [-wlx*0.5] [-wly*0.5] ...                        [-wlx*0.5] [-wly*0.5] [wlx*0.5] [-wly*0.5] ...                        [wlx*0.5] [-wly*0.5] [wlx*0.5] [-wly*0.5*1.5] id 2 



def add_bianjie    y_pos=-wly*0.5*1.2    loop while y_pos<wly*0.5*1.2         command            ball create position [-wlx*0.5-bianjie_rad] [y_pos] radius [bianjie_rad] group left_bianjie            ball create position [wlx*0.5 bianjie_rad] [y_pos] radius [bianjie_rad] group right_bianjie        endcommand        y_pos =2*bianjie_rad    endloopend@add_bianjie

contact groupbehavior andcmat default type ball-ball model rrlinear method deformability emod @emod_shiyang kratio 1.5 property fric 5.0 rr_fric 5.0 cmat add 2 model linearcbond  method deform emod 7e6 kratio 1.5 ...                                        property cb_tenf 1e300 cb_shearf 1e300 rgap [bianjie_rad*0.01] ...                                        range group left_bianjie cmat add 3 model linearcbond method deform emod 7e6 kratio 1.5 ...                                        property cb_tenf 1e300 cb_shearf 1e300 rgap [bianjie_rad*0.01] ...                                        range group right_bianjie cmat apply range group left_bianjiecmat apply range group right_bianjieball attribute density 1.5e3 damp 0.7 range group left_bianjieball attribute density 1.5e3 damp 0.7 range group right_bianjiecleancontact method bond gap [bianjie_rad*0.1][weiya=math.abs(txx)]def upadat_bianjie_force    loop foreach bp ball.groupmap("left_bianjie")        if ball.contactnum(bp)>=2 then            f=vector(0,0)            loop foreach ct ball.contactmap(bp)                if contact.model(ct)="linearcbond" then                    bp2=contact.end1(ct)                    if bp2=bp then                        bp2=contact.end2(ct)                    endif                    dirt=ball.pos(bp2)-ball.pos(bp)                    disp=math.sqrt(comp.x(dirt)^2 comp.y(dirt)^2)                    fmag=weiya*disp/2.0                     factor=1                    if ball.id(bp2)>ball.id(bp) then                        factor=-1                    endif                     norm=factor*vector(comp.y(dirt)*(-1),comp.x(dirt))/disp                     f= f norm*fmag                endif             endloop            ball.force.app(bp)=f        endif

    endloop     loop foreach bp ball.groupmap("right_bianjie")        if ball.contactnum(bp)>=2 then            f=vector(0,0)            loop foreach ct ball.contactmap(bp)                if contact.model(ct)="linearcbond" then                    bp2=contact.end1(ct)                    if bp2=bp then                        bp2=contact.end2(ct)                    endif                    dirt=ball.pos(bp2)-ball.pos(bp)                    disp=math.sqrt(comp.x(dirt)^2 comp.y(dirt)^2)                    fmag=weiya*disp/2.0                    factor=-1                    if ball.id(bp2)>ball.id(bp) then                        factor=1                    endif                    norm=factor*vector(comp.y(dirt)*(-1),comp.x(dirt))/disp                      f= f norm*fmag                endif            endloop            ball.force.app(bp)=f        endif

    endloopend

def initUpDownMo     loop foreach bp ball.groupmap("left_bianjie")        if ball.pos.y(bp)>wly*0.5-bianjie_rad then            ball.group(bp,10)="upMo"         endif        if ball.pos.y(bp)<-wly*0.5 bianjie_rad then            ball.group(bp,10)="downMo"         endif     endloop     loop foreach bp ball.groupmap("right_bianjie")        if ball.pos.y(bp)>wly*0.5-bianjie_rad then            ball.group(bp,10)="upMo"        endif        if ball.pos.y(bp)<-wly*0.5 bianjie_rad then            ball.group(bp,10)="downMo"                    endif     endloopend@initUpDownMo

ball fix vel spin range group upMo slot 10ball fix vel spin range group downMo slot 10



[now_timestep=global.step-1]def servo_bianjie    stress_jiance    if global.step > now_timestep then        upadat_bianjie_force        now_timestep=global.step freq    endifend

[wpDown=wall.find(2)][wpUp=wall.find(1)]

def sevro_walls    compute_stress   if timestepNow<global.step then        get_g(sevro_factor)        timestepNow =sevro_freq    endif        if do_ySevro=true then        Yvel=gy*(wyss-tyy)        wall.vel.y(wpUp)=-Yvel        wall.vel.y(wpDown)=Yvel        loop foreach bp ball.groupmap("upMo",10)            ball.vel.y(bp)=-Yvel        endloop        loop foreach bp ball.groupmap("downMo",10)            ball.vel.y(bp)=Yvel        endloop    endifend

def computer_chiCun    wly=wall.pos.y(wpUp)-wall.pos.y(wpDown)enddef compute_stress    computer_chiCun    wyss=-(wall.force.contact.y(wpUp)-wall.force.contact.y(wpDown))*0.5/wlxend

def get_g(fac)    gy=0    zongKNY=100e6*2*10    loop foreach ct wall.contactmap(wpUp)        zongKNY =contact.prop(ct,"kn")    endloop    loop foreach ct wall.contactmap(wpDown)        zongKNY =contact.prop(ct,"kn")    endloop    gy=fac*wlx/(zongKNY*global.timestep)    endset fish callback -1.0 @servo_bianjieset fish callback -1.0 @sevro_wallsmeasure create position 0 0 radius [wlx*0.3] id 1[mp=measure.find(1)]def stress_jiance    stressXX=measure.stress.xx(mp)end

history deletehistory id 1 @stressXXhistory id 2 @wysscycle 1solve

save rouxingmo

伺服结束后如图：

640 (1).png

可以看到是符合我们需求的，我们看一下颗粒受力图：

640 (2).png

当然我们还布置了测量圆测了一下横向应力：

横向应力距离目标应力有误差，来源于颗粒效应以及测量圆的误差。

三、加载

最后一步没啥新鲜的：

restore rouxingmo

wall attribute displacement multiply 0ball attribute displacement multiply 0set fish callback -1.0 remove @sevro_walls

[streainRatio=1e-2]

def get_ding_pos    ftdown=wall.facet.find(5)    ftup=wall.facet.find(2)    y_ding_pos=math.abs(wall.facet.pos.y(ftup))    end@get_ding_pos

wall attribute yvelocity [y_ding_pos*streainRatio] range id 2wall attribute yvelocity [-y_ding_pos*streainRatio] range id 1

ball attribute yvelocity [-y_ding_pos*streainRatio] range group upMo slot 10ball attribute yvelocity [y_ding_pos*streainRatio] range group downMo slot 10





[Iy0=y_ding_pos*2]def computer_stress_strain        q=wyss-txx     weyy=(-1)*(Iy0-(wall.facet.pos.y(ftup)-wall.facet.pos.y(ftdown)))/Iy0     wyss=-(wall.force.contact.y(wpUp)-wall.force.contact.y(wpDown))*0.5/wlxend

set fish callback -1.1 @computer_stress_strain

history deletehistory id 1 @wysshistory id 2 @weyyhistory id 3 @q

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

[baocunpinlv=2e-3][time_record=weyy 1][count=0]def savefile        if time_record-weyy >= 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