0 引言
这里想要尝试做的是同济大学黄雨老师的 Computational assessment of baffle performanceagainst rapid granular flows。 这篇文章从模拟方面论证了不同挡板结构对碎屑物冲击的支护性能评价。
1、成样
这里模拟颗粒运移,所以对初始的模型平衡要求并不是特别高,为了计算省事,只solve到-2次方就停止了,有计算条件的同学可以到-5次方。
我们的模型大概是下面这个样子:
这里的建模方式是在水平向进行墙体和颗粒区域的构建,然后将其旋转到指定角度。在初始模型中,颗粒的前端应当是有挡土结构的。
model new
def par
angle=40
rdmin=0.004
rdmax=rdmin*1.5
shiyang_width=0.44
shiyang_height=0.15
end
@par
geometry set "shiyang"
geometry generate box [-2.19-shiyang_width] -2.19 0 0.4 0 @shiyang_height
geometry rotate angle @angle axis (0,1,0) origin (0,0,0)
model domain extent [(-2.19-shiyang_width)*math.cos(angle*math.pi/180.0)] ...
[(-2.19)*math.cos(angle*math.pi/180.0) shiyang_height] 0 0.4 [(2.19)*math.sin(angle*math.pi/180.0)-shiyang_height] ...
[(2.19 shiyang_width)*math.sin(angle*math.pi/180.0) shiyang_height] ...
ball distribute radius @rdmin @rdmax porosity 0.4 range geometry-space "shiyang" ...
inside
model domain extent -4 1 -1 1 0 4
wall generate id 1 box -3 0 0 0.4 one-wall
wall generate id 2 polygon (-2.19,0,0) (-2.19,0.4,0) (-2.19,0.4,0.4) ...
(-2.19,0,0.4)
wall rotate axis (0,1,0) point (0,0,0) angle @angle
ball attribute density 2.7e3 damp 0.7
contact cmat default type ball-facet model linear method deform emod 100e6 kratio 1.5 property fric 0
contact cmat default type ball-ball model linear method deform emod 100e6 kratio 1.5 property fric 0.2
model cycle 2000 calm 50
ball delete range geometry-space "shiyang" outside
model gravity 9.8
model cycle 1
model solve ratio-average 1e-2
ball attribute damp 0.2
model cycle 1
model solve ratio-average 1e-2
model save "sample"
生成的模型如图:
可以看到在重力影响下,式样是不能维持一开始的方形。
2 下落
这个事情不难,只需要将初始的挡土结构删掉,式样便可以在重力作用下实现下落。需要注意的是,由于一开始我们的方形墙体指定的是onewall,所以无法从wall上面直接拿到冲击力,需要遍历面片上的接触来获得需要的力。这里是每0.5s保存一个sav文件,共计算了10s。
model restore "sample"
model mechanical time-total 0
def wall_init
wp5=wall.facet.find(5)
wp6=wall.facet.find(6)
end
@wall_init
def jiance
whilestepping
force=0
loop foreach ct wall.facet.contactmap(wp5)
force =math.mag(contact.force.global(ct))
endloop
loop foreach ct wall.facet.contactmap(wp6)
force =math.mag(contact.force.global(ct))
endloop
end
wall delete walls range id 2
[baocunpinlv=0.5]
[time_record=-10]
[count=0]
def savefile
time=mech.time.total
if time-time_record >= baocunpinlv then
filename=string.build("jieguo%1",count)
command
model save @filename
endcommand
time_record=time
count =1
endif
end
fish callback add @savefile -1.0
history id 1 @force
model cycle 1
model solve time 10
model save "result"
3 挡板结构模拟
挡板结构我是写了一个dat,然后在下落的时候call一下,便可以自动生产挡板了。这里比较的为方形和三角形挡板。
即:
挡板的逻辑和前面类似,也是先在x-y面上定义,然后旋转到指定角度上。
方形挡板:
[d=0.4/10.0]
[DBack=0.4/3.0]
geometry set "pile_back"
geometry generate box [-d*0.5] [d*0.5] [-d*0.5] [d*0.5] 0 0.3
geometry translate -0.42 0 0
geometry rotate angle @angle axis (0,1,0) origin (0,0,0)
wall import from-geometry "pile_back" name "pile_back_1"
def addpile_back
wall_count=2
loop n(1,3)
wall_name=string.build("pile_back_%1",wall_count)
command
geometry translate 0 @DBack 0
wall import from-geometry "pile_back" name @wall_name
endcommand
wall_count =1
endloop
end
@addpile_back
geometry set "pile_front"
geometry generate box [-d*0.5] [d*0.5] [-d*0.5] [d*0.5] 0 0.3
geometry translate -0.21 [DBack*0.5] 0
geometry rotate angle @angle axis (0,1,0) origin (0,0,0)
wall import from-geometry "pile_front" name "pile_front_1"
def addpile_front
wall_count=2
loop n(1,2)
wall_name=string.build("pile_front_%1",wall_count)
command
geometry translate 0 @DBack 0
wall import from-geometry "pile_front" name @wall_name
endcommand
wall_count =1
endloop
end
@addpile_front
三角形挡板:
[d=0.4/10.0]
[DBack=0.4/3.0]
geometry set "pile_back"
geometry edge create by-position (0,[-d*0.5],0) ([-d*0.5],0,0) (0,[d*0.5],0) (0,[-d*0.5],0)
geometry generate from-edges extrude (0,0,0.3) segments 1
geometry triangulate clean
geometry translate -0.42 0 0
geometry rotate angle @angle axis (0,1,0) origin (0,0,0)
wall import from-geometry "pile_back" name "pile_back_1"
def addpile_back
wall_count=2
loop n(1,3)
wall_name=string.build("pile_back_%1",wall_count)
command
geometry translate 0 @DBack 0
wall import from-geometry "pile_back" name @wall_name
endcommand
wall_count =1
endloop
end
@addpile_back
geometry set "pile_front"
geometry edge create by-position (0,[-d*0.5],0) ([-d*0.5],0,0) (0,[d*0.5],0) (0,[-d*0.5],0)
geometry generate from-edges extrude (0,0,0.3) segments 1
geometry triangulate clean
geometry translate -0.21 [DBack*0.5] 0
geometry rotate angle @angle axis (0,1,0) origin (0,0,0)
wall import from-geometry "pile_front" name "pile_front_1"
def addpile_front
wall_count=2
loop n(1,2)
wall_name=string.build("pile_front_%1",wall_count)
command
geometry translate 0 @DBack 0
wall import from-geometry "pile_front" name @wall_name
endcommand
wall_count =1
endloop
end
@addpile_front
模型为:
4 结果分析
这种颗粒运移的图肯定是以动图为主了。
无挡板
侧面:
正面:
方形挡板
侧面:
正面:
三角形挡板
侧面:
正面:
通过上面的动图可以非常明显的感觉到挡板结构对后一段的碎屑物有非常突出的支护性能。
这里看一下运行到5s的时候的力链图:
无挡板:
这里可以明显的看出,前缘的碎屑物和后缘有力学上的不连续性。
方形挡板:
三角形挡板:
可以看到挡板间形成的土拱有效的抵抗了后缘碎屑物的冲击。
这里再看一下冲击力的大小,下图为三角形挡板支护时的冲击力,由于初始的颗粒没有阻碍,所以初期会有很多颗粒以比较大的速度冲击顶板上,导致初期有比较大的波动。这里用FFT每50点进行平滑处理,处理后的为红色线。
因为数据量比较大,origin没抗住,这里用excel简单画了一下:
可以看出挡板结构对碎屑流有比较好的支护性能,方形挡板略优于三角形,但是三角形耗材少,具体工程指导得按实际来。