最近旁边的小伙伴在做有限元的动力问题,需要用到粘弹性边界,于是我也使用离散元模拟了一下。这里采用三种边界进行对比:固定边界、墙体阻尼边界、自由阻尼边界。
一、模型工况
首先成一个矩形的样,长宽为1.5*0.75,四周墙体速度固定,进行平衡。代码如下:
new
set random 10001
def par
width=1.5
height=width*0.5
poro=0.08
emod_shiyang=100e6
end
@par
domain extent [-height*2] [height*2]
wall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] expand 1.5
cmat default model linear method deformability emod @emod_shiyang kratio 1.5
ball 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.7
cycle 1000 calm 10
solve
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 sample
成样后如下图所示:
之后给左边墙体加一个激振,中间布置三个测点,测试其时程曲线。
入射边界先振动一个周期,后面停止振动,运行五个周期,于是一个周期后的测点振动都是反射引起的。时程图如下:
二、边界处理
1、固定边界
这个很简单,什么都不用做,直接振动即可:
restore sample
ball property fric 0.5
ball attribute damp 0
set mech age 0
def addcedian
cedian1=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
@addcedian
def get_vel
whilestepping
sudu1=ball.vel.x(cedian1)
sudu2=ball.vel.x(cedian2)
sudu3=ball.vel.x(cedian3)
time=mech.age
end
@get_vel
[wp=wall.find(4)]
[fengzhi=width*0.2]
[pinlv=20]
def setvel
vel=fengzhi*math.sin(2*math.pi*pinlv*time)
wall.vel.x(wp)=vel
end
@setvel
history id 1 @vel
history id 2 @sudu1
history id 3 @sudu2
history id 4 @sudu3
history id 5 @time
set fish callback -1.0 @setvel
solve time [1/float(pinlv)]
set fish callback -1.0 remove @setvel
wall attribute vel 0
solve time [5/float(pinlv)]
save nodamp
2、墙体阻尼边界
这个边界处理只需要给ball-facet接触加阻尼系数就可以了:
restore sample
ball property fric 0.5
ball attribute damp 0
set mech age 0
cmat 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 apply
def addcedian
cedian1=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
@addcedian
def get_vel
whilestepping
time=mech.age
sudu1=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 setvel
vel=fengzhi*math.sin(2*math.pi*pinlv*time)
wall.vel.x(wp)=vel
end
@setvel
history id 1 @vel
history id 2 @sudu1
history id 3 @sudu2
history id 4 @sudu3
history id 5 @time
set fish callback -1.0 @setvel
solve time [1/float(pinlv)]
set fish callback -1.0 remove @setvel
wall attribute vel 0
solve time [5/float(pinlv)]
save walldamp
3、自由阻尼边界
这个说法有很多,也有说是人工阻尼边界,说法无所谓,关键是原理。这部分原理是加一个可以自由移动的边界,这个边界要同时满足两个条件:1)满足平衡条件;2)阻碍质点位移。
从这两点出发,在PFC中施加自由阻尼边界分两步走,首先将颗粒与墙的力施加在边界颗粒上,同时删除墙;之后根据颗粒为速度施加一定的阻尼力。
对于岩石来说,这个比较容易一点,只需要给边界颗粒加反力,然后设置边界颗粒的阻尼系数即可。但是对于砂土而言便困难一点,因为砂土有时候会从边界颗粒缝隙挤出。
于是我在这里提出一个关于砂土自由阻尼边界的做法,使用clump来施加反力以及阻尼。做法是将墙的接触力给clump,同时删除墙,给clump一定的阻尼系数。
下面为阻尼边界的施加:
restore sample
[jiazaibanhoudu=width*0.01]
geometry set shujiazaiban
geometry generate box [-jiazaibanhoudu] [0] ...
[-height*0.5*1.5] [height*0.5*1.5]
geometry set hengjiazaiban
geometry 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 ...
surfcalculate
clump 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 ...
surfcalculate
clump replicate id 2 name shujiazaiban ...
x [width*0.5 jiazaibanhoudu*0.5] y 0
def addforce
wpright=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
@addforce
cmat add 1 model null range contact type pebble-facet
cmat add 2 model null range contact type pebble-pebble
clump attribute damp 0.7
clump fix yvel spin range id 2
clump fix xvel spin range id 1 3
wall delete walls range id 1
wall delete walls range id 2
wall delete walls range id 3
cycle 1
solve
save addbianjie
之后是加载测量:
restore addbianjie
contact property dp_nratio 0.7 dp_sratio 0.7 range contact type ball-pebble
clump attribute density 1e3
ball property fric 0.5
ball attribute damp 0
set mech age 0
def addcedian
cedian1=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
@addcedian
def get_vel
whilestepping
time=mech.age
sudu1=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 setvel
vel=fengzhi*math.sin(2*math.pi*pinlv*time)
wall.vel.x(wp)=vel
end
@setvel
history id 1 @vel
history id 2 @sudu1
history id 3 @sudu2
history id 4 @sudu3
history id 5 @time
set fish callback -1.0 @setvel
solve time [1/float(pinlv)]
set fish callback -1.0 remove @setvel
wall attribute vel 0
solve time [5/float(pinlv)]
save freedamp
三、结果对比
下面对比一下三个方法的结果。
1、固定边界
可以看出一次振动后,内部颗粒做等幅的振动,没有衰减现象产生。
2、墙体阻尼边界
相比于固定边界,振动有了明显的 ,说明采用墙体阻尼边界处理动力问题是可行的。
3、自由阻尼边界
效果来说,其和墙体阻尼差不多。
四、总结
粘弹性边界是有限元方法中处理动力边界的方法,这是有限元处理动力问题的缺陷算是,也就是节点作为边界。但是在离散元中,边界是由wall来充当,并且wall和颗粒可以分开,也可以设置各种弹性和粘性参数。从粘弹性边界的原理出发,墙体阻尼边界足够满足粘弹性边界的两个条件了。所以这里建议采用墙体边界便足够进行动力问题分析了。