0 引言
这里想要尝试做的是同济大学黄雨老师的 Computational assessment of baffle performanceagainst rapid granular flows。 这篇文章从模拟方面论证了不同挡板结构对碎屑物冲击的支护性能评价。
1、成样
这里模拟颗粒运移,所以对初始的模型平衡要求并不是特别高,为了计算省事,只solve到-2次方就停止了,有计算条件的同学可以到-5次方。
我们的模型大概是下面这个样子:

这里的建模方式是在水平向进行墙体和颗粒区域的构建,然后将其旋转到指定角度。在初始模型中,颗粒的前端应当是有挡土结构的。
model newdef parangle=40rdmin=0.004rdmax=rdmin*1.5shiyang_width=0.44shiyang_height=0.15end@pargeometry set "shiyang"geometry generate box [-2.19-shiyang_width] -2.19 0 0.4 0 @shiyang_heightgeometry 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" ...insidemodel domain extent -4 1 -1 1 0 4wall generate id 1 box -3 0 0 0.4 one-wallwall 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 @angleball attribute density 2.7e3 damp 0.7contact cmat default type ball-facet model linear method deform emod 100e6 kratio 1.5 property fric 0contact cmat default type ball-ball model linear method deform emod 100e6 kratio 1.5 property fric 0.2model cycle 2000 calm 50ball delete range geometry-space "shiyang" outsidemodel gravity 9.8model cycle 1model solve ratio-average 1e-2ball attribute damp 0.2model cycle 1model solve ratio-average 1e-2model save "sample"
生成的模型如图:

可以看到在重力影响下,式样是不能维持一开始的方形。
2 下落
这个事情不难,只需要将初始的挡土结构删掉,式样便可以在重力作用下实现下落。需要注意的是,由于一开始我们的方形墙体指定的是onewall,所以无法从wall上面直接拿到冲击力,需要遍历面片上的接触来获得需要的力。这里是每0.5s保存一个sav文件,共计算了10s。
model restore "sample"model mechanical time-total 0def wall_initwp5=wall.facet.find(5)wp6=wall.facet.find(6)end@wall_initdef jiancewhilesteppingforce=0loop foreach ct wall.facet.contactmap(wp5)force =math.mag(contact.force.global(ct))endlooploop foreach ct wall.facet.contactmap(wp6)force =math.mag(contact.force.global(ct))endloopendwall delete walls range id 2[baocunpinlv=0.5][time_record=-10][count=0]def savefiletime=mech.time.totalif time-time_record >= baocunpinlv thenfilename=string.build("jieguo%1",count)commandmodel save @filenameendcommandtime_record=timecount =1endifendfish callback add @savefile -1.0history id 1 @forcemodel cycle 1model solve time 10model 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.3geometry translate -0.42 0 0geometry rotate angle @angle axis (0,1,0) origin (0,0,0)wall import from-geometry "pile_back" name "pile_back_1"def addpile_backwall_count=2loop n(1,3)wall_name=string.build("pile_back_%1",wall_count)commandgeometry translate 0 @DBack 0wall import from-geometry "pile_back" name @wall_nameendcommandwall_count =1endloopend@addpile_backgeometry set "pile_front"geometry generate box [-d*0.5] [d*0.5] [-d*0.5] [d*0.5] 0 0.3geometry translate -0.21 [DBack*0.5] 0geometry rotate angle @angle axis (0,1,0) origin (0,0,0)wall import from-geometry "pile_front" name "pile_front_1"def addpile_frontwall_count=2loop n(1,2)wall_name=string.build("pile_front_%1",wall_count)commandgeometry translate 0 @DBack 0wall import from-geometry "pile_front" name @wall_nameendcommandwall_count =1endloopend@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 1geometry triangulate cleangeometry translate -0.42 0 0geometry rotate angle @angle axis (0,1,0) origin (0,0,0)wall import from-geometry "pile_back" name "pile_back_1"def addpile_backwall_count=2loop n(1,3)wall_name=string.build("pile_back_%1",wall_count)commandgeometry translate 0 @DBack 0wall import from-geometry "pile_back" name @wall_nameendcommandwall_count =1endloopend@addpile_backgeometry 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 1geometry triangulate cleangeometry translate -0.21 [DBack*0.5] 0geometry rotate angle @angle axis (0,1,0) origin (0,0,0)wall import from-geometry "pile_front" name "pile_front_1"def addpile_frontwall_count=2loop n(1,2)wall_name=string.build("pile_front_%1",wall_count)commandgeometry translate 0 @DBack 0wall import from-geometry "pile_front" name @wall_nameendcommandwall_count =1endloopend@addpile_front
模型为:

4 结果分析
这种颗粒运移的图肯定是以动图为主了。
无挡板
侧面:
.gif?imageView2/1/w/631/h/544)
正面:
.gif?imageView2/1/w/718/h/534)
方形挡板
侧面:
.gif?imageView2/1/w/653/h/566)
正面:
.gif?imageView2/1/w/725/h/548)
三角形挡板
侧面:
.gif?imageView2/1/w/733/h/523)
正面:
.gif?imageView2/1/w/777/h/589)
通过上面的动图可以非常明显的感觉到挡板结构对后一段的碎屑物有非常突出的支护性能。
这里看一下运行到5s的时候的力链图:
无挡板:

这里可以明显的看出,前缘的碎屑物和后缘有力学上的不连续性。
方形挡板:

三角形挡板:

可以看到挡板间形成的土拱有效的抵抗了后缘碎屑物的冲击。
这里再看一下冲击力的大小,下图为三角形挡板支护时的冲击力,由于初始的颗粒没有阻碍,所以初期会有很多颗粒以比较大的速度冲击顶板上,导致初期有比较大的波动。这里用FFT每50点进行平滑处理,处理后的为红色线。

因为数据量比较大,origin没抗住,这里用excel简单画了一下:

可以看出挡板结构对碎屑流有比较好的支护性能,方形挡板略优于三角形,但是三角形耗材少,具体工程指导得按实际来。