最近旁边的小伙伴在做有限元的动力问题,需要用到粘弹性边界,于是我也使用离散元模拟了一下。这里采用三种边界进行对比:固定边界、墙体阻尼边界、自由阻尼边界。
一、模型工况
首先成一个矩形的样,长宽为1.5*0.75,四周墙体速度固定,进行平衡。代码如下:
newset random 10001def parwidth=1.5height=width*0.5poro=0.08emod_shiyang=100e6end@pardomain extent [-height*2] [height*2]wall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] expand 1.5cmat default model linear method deformability emod @emod_shiyang kratio 1.5ball distribute porosity @poro radius 0.006 0.009 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 10solveball 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 sample
成样后如下图所示:

之后给左边墙体加一个激振,中间布置三个测点,测试其时程曲线。

入射边界先振动一个周期,后面停止振动,运行五个周期,于是一个周期后的测点振动都是反射引起的。时程图如下:

二、边界处理
1、固定边界
这个很简单,什么都不用做,直接振动即可:
restore sampleball property fric 0.5ball attribute damp 0set mech age 0def addcediancedian1=ball.near(-width*0.4,0)cedian2=ball.near(0,0)cedian3=ball.near(width*0.4,0)ball.group(cedian1)="cedian"ball.group(cedian2)="cedian"ball.group(cedian3)="cedian"end@addcediandef get_velwhilesteppingsudu1=ball.vel.x(cedian1)sudu2=ball.vel.x(cedian2)sudu3=ball.vel.x(cedian3)time=mech.ageend@get_vel[wp=wall.find(4)][fengzhi=width*0.2][pinlv=20]def setvelvel=fengzhi*math.sin(2*math.pi*pinlv*time)wall.vel.x(wp)=velend@setvelhistory id 1 @velhistory id 2 @sudu1history id 3 @sudu2history id 4 @sudu3history id 5 @timeset fish callback -1.0 @setvelsolve time [1/float(pinlv)]set fish callback -1.0 remove @setvelwall attribute vel 0solve time [5/float(pinlv)]save nodamp
2、墙体阻尼边界
这个边界处理只需要给ball-facet接触加阻尼系数就可以了:
restore sampleball property fric 0.5ball attribute damp 0set mech age 0cmat add 1 model linear method deformability emod @emod_shiyang kratio 1.5 ...property dp_nratio 2.0 dp_sratio 2.0 range contact type ball-facet x [-width*0.4] [width]cmat applydef addcediancedian1=ball.near(-width*0.4,0)cedian2=ball.near(0,0)cedian3=ball.near(width*0.4,0)ball.group(cedian1)="cedian"ball.group(cedian2)="cedian"ball.group(cedian3)="cedian"end@addcediandef get_velwhilesteppingtime=mech.agesudu1=ball.vel.x(cedian1)sudu2=ball.vel.x(cedian2)sudu3=ball.vel.x(cedian3)end@get_vel[wp=wall.find(4)][fengzhi=width*0.3][pinlv=10]def setvelvel=fengzhi*math.sin(2*math.pi*pinlv*time)wall.vel.x(wp)=velend@setvelhistory id 1 @velhistory id 2 @sudu1history id 3 @sudu2history id 4 @sudu3history id 5 @timeset fish callback -1.0 @setvelsolve time [1/float(pinlv)]set fish callback -1.0 remove @setvelwall attribute vel 0solve time [5/float(pinlv)]save walldamp
3、自由阻尼边界
这个说法有很多,也有说是人工阻尼边界,说法无所谓,关键是原理。这部分原理是加一个可以自由移动的边界,这个边界要同时满足两个条件:1)满足平衡条件;2)阻碍质点位移。
从这两点出发,在PFC中施加自由阻尼边界分两步走,首先将颗粒与墙的力施加在边界颗粒上,同时删除墙;之后根据颗粒为速度施加一定的阻尼力。
对于岩石来说,这个比较容易一点,只需要给边界颗粒加反力,然后设置边界颗粒的阻尼系数即可。但是对于砂土而言便困难一点,因为砂土有时候会从边界颗粒缝隙挤出。
于是我在这里提出一个关于砂土自由阻尼边界的做法,使用clump来施加反力以及阻尼。做法是将墙的接触力给clump,同时删除墙,给clump一定的阻尼系数。
下面为阻尼边界的施加:
restore sample[jiazaibanhoudu=width*0.01]geometry set shujiazaibangeometry generate box [-jiazaibanhoudu] [0] ...[-height*0.5*1.5] [height*0.5*1.5]geometry set hengjiazaibangeometry generate box [-width*0.5*1.5] [width*0.5*1.5] [-jiazaibanhoudu] [0]clump template create name hengjiazaiban ...geometry hengjiazaiban ...bubblepack ratio 0.2 distance 160 ...surfcalculateclump replicate id 1 name hengjiazaiban ...x 0 y [-height*0.5-jiazaibanhoudu*0.5]clump replicate id 3 name hengjiazaiban ...x 0 y [height*0.5 jiazaibanhoudu*0.5]clump template create name shujiazaiban ...geometry shujiazaiban ...bubblepack ratio 0.2 distance 160 ...surfcalculateclump replicate id 2 name shujiazaiban ...x [width*0.5 jiazaibanhoudu*0.5] y 0def addforcewpright=wall.find(2)cpright=clump.find(2)clump.force.app.x(cpright)=(-1)*wall.force.contact.x(wpright)wpup=wall.find(3)cpup=clump.find(3)clump.force.app.y(cpup)=(-1)*wall.force.contact.y(wpup)wpdown=wall.find(1)cpdown=clump.find(1)clump.force.app.y(cpdown)=(-1)*wall.force.contact.y(wpdown)end@addforcecmat add 1 model null range contact type pebble-facetcmat add 2 model null range contact type pebble-pebbleclump attribute damp 0.7clump fix yvel spin range id 2clump fix xvel spin range id 1 3wall delete walls range id 1wall delete walls range id 2wall delete walls range id 3cycle 1solvesave addbianjie
之后是加载测量:
restore addbianjiecontact property dp_nratio 0.7 dp_sratio 0.7 range contact type ball-pebbleclump attribute density 1e3ball property fric 0.5ball attribute damp 0set mech age 0def addcediancedian1=ball.near(-width*0.4,0)cedian2=ball.near(0,0)cedian3=ball.near(width*0.4,0)ball.group(cedian1)="cedian"ball.group(cedian2)="cedian"ball.group(cedian3)="cedian"end@addcediandef get_velwhilesteppingtime=mech.agesudu1=ball.vel.x(cedian1)sudu2=ball.vel.x(cedian2)sudu3=ball.vel.x(cedian3)end@get_vel[wp=wall.find(4)][fengzhi=width*0.2][pinlv=20]def setvelvel=fengzhi*math.sin(2*math.pi*pinlv*time)wall.vel.x(wp)=velend@setvelhistory id 1 @velhistory id 2 @sudu1history id 3 @sudu2history id 4 @sudu3history id 5 @timeset fish callback -1.0 @setvelsolve time [1/float(pinlv)]set fish callback -1.0 remove @setvelwall attribute vel 0solve time [5/float(pinlv)]save freedamp
三、结果对比
下面对比一下三个方法的结果。
1、固定边界
可以看出一次振动后,内部颗粒做等幅的振动,没有衰减现象产生。

2、墙体阻尼边界
相比于固定边界,振动有了明显的 ,说明采用墙体阻尼边界处理动力问题是可行的。

3、自由阻尼边界
效果来说,其和墙体阻尼差不多。

四、总结
粘弹性边界是有限元方法中处理动力边界的方法,这是有限元处理动力问题的缺陷算是,也就是节点作为边界。但是在离散元中,边界是由wall来充当,并且wall和颗粒可以分开,也可以设置各种弹性和粘性参数。从粘弹性边界的原理出发,墙体阻尼边界足够满足粘弹性边界的两个条件了。所以这里建议采用墙体边界便足够进行动力问题分析了。